1 |
|
|
! |
2 |
|
|
MODULE slab_heat_transp_mod |
3 |
|
|
! |
4 |
|
|
! Slab ocean : temperature tendencies due to horizontal diffusion |
5 |
|
|
! and / or Ekman transport |
6 |
|
|
|
7 |
|
|
USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo |
8 |
|
|
IMPLICIT NONE |
9 |
|
|
|
10 |
|
|
! Variables copied over from dyn3d dynamics: |
11 |
|
|
REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area |
12 |
|
|
!$OMP THREADPRIVATE(fext) |
13 |
|
|
REAL,SAVE,ALLOCATABLE :: beta(:) ! df/dy |
14 |
|
|
!$OMP THREADPRIVATE(beta) |
15 |
|
|
REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area |
16 |
|
|
!$OMP THREADPRIVATE(unsairez) |
17 |
|
|
REAL,SAVE,ALLOCATABLE :: unsaire(:) |
18 |
|
|
!$OMP THREADPRIVATE(unsaire) |
19 |
|
|
REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m) |
20 |
|
|
!$OMP THREADPRIVATE(cu) |
21 |
|
|
REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m) |
22 |
|
|
!$OMP THREADPRIVATE(cv) |
23 |
|
|
REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points) |
24 |
|
|
!$OMP THREADPRIVATE(cuvsurcv) |
25 |
|
|
REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points) |
26 |
|
|
!$OMP THREADPRIVATE(cvusurcu) |
27 |
|
|
REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area |
28 |
|
|
!$OMP THREADPRIVATE(aire) |
29 |
|
|
REAL,SAVE :: apoln ! area of north pole points |
30 |
|
|
!$OMP THREADPRIVATE(apoln) |
31 |
|
|
REAL,SAVE :: apols ! area of south pole points |
32 |
|
|
!$OMP THREADPRIVATE(apols) |
33 |
|
|
REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells |
34 |
|
|
!$OMP THREADPRIVATE(aireu) |
35 |
|
|
REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells |
36 |
|
|
!$OMP THREADPRIVATE(airev) |
37 |
|
|
|
38 |
|
|
! Local parameters for slab transport |
39 |
|
|
LOGICAL,SAVE :: alpha_var ! variable coef for deep temp (1 layer) |
40 |
|
|
!$OMP THREADPRIVATE(alpha_var) |
41 |
|
|
LOGICAL,SAVE :: slab_upstream ! upstream scheme ? (1 layer) |
42 |
|
|
!$OMP THREADPRIVATE(slab_upstream) |
43 |
|
|
LOGICAL,SAVE :: slab_sverdrup ! use wind stress curl at equator |
44 |
|
|
!$OMP THREADPRIVATE(slab_sverdrup) |
45 |
|
|
LOGICAL,SAVE :: slab_tyeq ! use merid wind stress at equator |
46 |
|
|
!$OMP THREADPRIVATE(slab_tyeq) |
47 |
|
|
LOGICAL,SAVE :: ekman_zonadv ! use zonal advection by Ekman currents |
48 |
|
|
!$OMP THREADPRIVATE(ekman_zonadv) |
49 |
|
|
LOGICAL,SAVE :: ekman_zonavg ! zonally average wind stress |
50 |
|
|
!$OMP THREADPRIVATE(ekman_zonavg) |
51 |
|
|
|
52 |
|
|
REAL,SAVE :: alpham |
53 |
|
|
!$OMP THREADPRIVATE(alpham) |
54 |
|
|
REAL,SAVE :: gmkappa |
55 |
|
|
!$OMP THREADPRIVATE(gmkappa) |
56 |
|
|
REAL,SAVE :: gm_smax |
57 |
|
|
!$OMP THREADPRIVATE(gm_smax) |
58 |
|
|
|
59 |
|
|
! geometry variables : f, beta, mask... |
60 |
|
|
REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux |
61 |
|
|
!$OMP THREADPRIVATE(zmasqu) |
62 |
|
|
REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux |
63 |
|
|
!$OMP THREADPRIVATE(zmasqv) |
64 |
|
|
REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points |
65 |
|
|
!$OMP THREADPRIVATE(unsfv) |
66 |
|
|
REAL,SAVE,ALLOCATABLE :: unsbv(:) ! 1/beta |
67 |
|
|
!$OMP THREADPRIVATE(unsbv) |
68 |
|
|
REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag) |
69 |
|
|
!$OMP THREADPRIVATE(unsev) |
70 |
|
|
REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points |
71 |
|
|
!$OMP THREADPRIVATE(unsfu) |
72 |
|
|
REAL,SAVE,ALLOCATABLE :: unseu(:) |
73 |
|
|
!$OMP THREADPRIVATE(unseu) |
74 |
|
|
|
75 |
|
|
! Routines from dyn3d, valid on global dynamics grid only: |
76 |
|
|
PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid |
77 |
|
|
PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points |
78 |
|
|
PRIVATE :: grad,diverg |
79 |
|
|
|
80 |
|
|
CONTAINS |
81 |
|
|
|
82 |
|
|
SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,& |
83 |
|
|
cu_,cuvsurcv_,cv_,cvusurcu_, & |
84 |
|
|
aire_,apoln_,apols_, & |
85 |
|
|
aireu_,airev_,rlatv, rad, omeg) |
86 |
|
|
! number of points in lon, lat |
87 |
|
|
IMPLICIT NONE |
88 |
|
|
! Routine copies some geometry variables from the dynamical core |
89 |
|
|
! see global vars for meaning |
90 |
|
|
INTEGER,INTENT(IN) :: ip1jm |
91 |
|
|
INTEGER,INTENT(IN) :: ip1jmp1 |
92 |
|
|
REAL,INTENT(IN) :: unsairez_(ip1jm) |
93 |
|
|
REAL,INTENT(IN) :: fext_(ip1jm) |
94 |
|
|
REAL,INTENT(IN) :: unsaire_(ip1jmp1) |
95 |
|
|
REAL,INTENT(IN) :: cu_(ip1jmp1) |
96 |
|
|
REAL,INTENT(IN) :: cuvsurcv_(ip1jm) |
97 |
|
|
REAL,INTENT(IN) :: cv_(ip1jm) |
98 |
|
|
REAL,INTENT(IN) :: cvusurcu_(ip1jmp1) |
99 |
|
|
REAL,INTENT(IN) :: aire_(ip1jmp1) |
100 |
|
|
REAL,INTENT(IN) :: apoln_ |
101 |
|
|
REAL,INTENT(IN) :: apols_ |
102 |
|
|
REAL,INTENT(IN) :: aireu_(ip1jmp1) |
103 |
|
|
REAL,INTENT(IN) :: airev_(ip1jm) |
104 |
|
|
REAL,INTENT(IN) :: rlatv(nbp_lat-1) |
105 |
|
|
REAL,INTENT(IN) :: rad |
106 |
|
|
REAL,INTENT(IN) :: omeg |
107 |
|
|
|
108 |
|
|
CHARACTER (len = 20) :: modname = 'slab_heat_transp' |
109 |
|
|
CHARACTER (len = 80) :: abort_message |
110 |
|
|
|
111 |
|
|
! Sanity check on dimensions |
112 |
|
|
if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. & |
113 |
|
|
(ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then |
114 |
|
|
abort_message="ini_slab_transp_geom Error: wrong array sizes" |
115 |
|
|
CALL abort_physic(modname,abort_message,1) |
116 |
|
|
endif |
117 |
|
|
! Allocations could be done only on master process/thread... |
118 |
|
|
allocate(unsairez(ip1jm)) |
119 |
|
|
unsairez(:)=unsairez_(:) |
120 |
|
|
allocate(fext(ip1jm)) |
121 |
|
|
fext(:)=fext_(:) |
122 |
|
|
allocate(unsaire(ip1jmp1)) |
123 |
|
|
unsaire(:)=unsaire_(:) |
124 |
|
|
allocate(cu(ip1jmp1)) |
125 |
|
|
cu(:)=cu_(:) |
126 |
|
|
allocate(cuvsurcv(ip1jm)) |
127 |
|
|
cuvsurcv(:)=cuvsurcv_(:) |
128 |
|
|
allocate(cv(ip1jm)) |
129 |
|
|
cv(:)=cv_(:) |
130 |
|
|
allocate(cvusurcu(ip1jmp1)) |
131 |
|
|
cvusurcu(:)=cvusurcu_(:) |
132 |
|
|
allocate(aire(ip1jmp1)) |
133 |
|
|
aire(:)=aire_(:) |
134 |
|
|
apoln=apoln_ |
135 |
|
|
apols=apols_ |
136 |
|
|
allocate(aireu(ip1jmp1)) |
137 |
|
|
aireu(:)=aireu_(:) |
138 |
|
|
allocate(airev(ip1jm)) |
139 |
|
|
airev(:)=airev_(:) |
140 |
|
|
allocate(beta(nbp_lat-1)) |
141 |
|
|
beta(:)=2*omeg*cos(rlatv(:))/rad |
142 |
|
|
|
143 |
|
|
END SUBROUTINE ini_slab_transp_geom |
144 |
|
|
|
145 |
|
|
SUBROUTINE ini_slab_transp(zmasq) |
146 |
|
|
|
147 |
|
|
! USE ioipsl_getin_p_mod, only: getin_p |
148 |
|
|
USE IOIPSL, ONLY : getin |
149 |
|
|
IMPLICIT NONE |
150 |
|
|
|
151 |
|
|
REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent |
152 |
|
|
REAL zmasq_2d((nbp_lon+1)*nbp_lat) |
153 |
|
|
REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter |
154 |
|
|
REAL eps ! epsilon friction timescale (s-1) |
155 |
|
|
INTEGER :: slab_ekman |
156 |
|
|
INTEGER i |
157 |
|
|
INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1 |
158 |
|
|
|
159 |
|
|
! Some definition for grid size |
160 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) |
161 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat |
162 |
|
|
iim=nbp_lon |
163 |
|
|
iip1=nbp_lon+1 |
164 |
|
|
jjp1=nbp_lat |
165 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) |
166 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat |
167 |
|
|
|
168 |
|
|
! Options for Heat transport |
169 |
|
|
! Alpha variable? |
170 |
|
|
alpha_var=.FALSE. |
171 |
|
|
CALL getin('slab_alphav',alpha_var) |
172 |
|
|
print *,'alpha variable',alpha_var |
173 |
|
|
! centered ou upstream scheme for meridional transport |
174 |
|
|
slab_upstream=.FALSE. |
175 |
|
|
CALL getin('slab_upstream',slab_upstream) |
176 |
|
|
print *,'upstream slab scheme',slab_upstream |
177 |
|
|
! Sverdrup balance at equator ? |
178 |
|
|
slab_sverdrup=.TRUE. |
179 |
|
|
CALL getin('slab_sverdrup',slab_sverdrup) |
180 |
|
|
print *,'Sverdrup balance',slab_sverdrup |
181 |
|
|
! Use tauy for meridional flux at equator ? |
182 |
|
|
slab_tyeq=.TRUE. |
183 |
|
|
CALL getin('slab_tyeq',slab_tyeq) |
184 |
|
|
print *,'Tauy forcing at equator',slab_tyeq |
185 |
|
|
! Use tauy for meridional flux at equator ? |
186 |
|
|
ekman_zonadv=.TRUE. |
187 |
|
|
CALL getin('slab_ekman_zonadv',ekman_zonadv) |
188 |
|
|
print *,'Use Ekman flow in zonal direction',ekman_zonadv |
189 |
|
|
! Use tauy for meridional flux at equator ? |
190 |
|
|
ekman_zonavg=.FALSE. |
191 |
|
|
CALL getin('slab_ekman_zonavg',ekman_zonavg) |
192 |
|
|
print *,'Use zonally-averaged wind stress ?',ekman_zonavg |
193 |
|
|
! Value of alpha |
194 |
|
|
alpham=2./3. |
195 |
|
|
CALL getin('slab_alpha',alpham) |
196 |
|
|
print *,'slab_alpha',alpham |
197 |
|
|
! GM k coefficient (m2/s) for 2-layers |
198 |
|
|
gmkappa=1000. |
199 |
|
|
CALL getin('slab_gmkappa',gmkappa) |
200 |
|
|
print *,'slab_gmkappa',gmkappa |
201 |
|
|
! GM k coefficient (m2/s) for 2-layers |
202 |
|
|
gm_smax=2e-3 |
203 |
|
|
CALL getin('slab_gm_smax',gm_smax) |
204 |
|
|
print *,'slab_gm_smax',gm_smax |
205 |
|
|
! ----------------------------------------------------------- |
206 |
|
|
! Define ocean / continent mask (no flux into continent cell) |
207 |
|
|
allocate(zmasqu(ip1jmp1)) |
208 |
|
|
allocate(zmasqv(ip1jm)) |
209 |
|
|
zmasqu=1. |
210 |
|
|
zmasqv=1. |
211 |
|
|
|
212 |
|
|
! convert mask to 2D grid |
213 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d) |
214 |
|
|
! put flux mask to 0 at boundaries of continent cells |
215 |
|
|
DO i=1,ip1jmp1-1 |
216 |
|
|
IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN |
217 |
|
|
zmasqu(i)=0. |
218 |
|
|
ENDIF |
219 |
|
|
END DO |
220 |
|
|
DO i=iip1,ip1jmp1,iip1 |
221 |
|
|
zmasqu(i)=zmasqu(i-iim) |
222 |
|
|
END DO |
223 |
|
|
DO i=1,ip1jm |
224 |
|
|
IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN |
225 |
|
|
zmasqv(i)=0. |
226 |
|
|
END IF |
227 |
|
|
END DO |
228 |
|
|
|
229 |
|
|
! ----------------------------------------------------------- |
230 |
|
|
! Coriolis and friction for Ekman transport |
231 |
|
|
slab_ekman=2 |
232 |
|
|
CALL getin("slab_ekman",slab_ekman) |
233 |
|
|
IF (slab_ekman.GT.0) THEN |
234 |
|
|
allocate(unsfv(ip1jm)) |
235 |
|
|
allocate(unsev(ip1jm)) |
236 |
|
|
allocate(unsfu(ip1jmp1)) |
237 |
|
|
allocate(unseu(ip1jmp1)) |
238 |
|
|
allocate(unsbv(ip1jm)) |
239 |
|
|
|
240 |
|
|
eps=1e-5 ! Drag |
241 |
|
|
CALL getin('slab_eps',eps) |
242 |
|
|
print *,'epsilon=',eps |
243 |
|
|
ff=fext*unsairez ! Coriolis |
244 |
|
|
! coefs to convert tau_x, tau_y to Ekman mass fluxes |
245 |
|
|
! on 2D grid v points |
246 |
|
|
! Compute correction factor [0 1] near the equator (f<<eps) |
247 |
|
|
IF (slab_sverdrup) THEN |
248 |
|
|
! New formulation, sharper near equator, when eps gives Rossby Radius |
249 |
|
|
DO i=1,ip1jm |
250 |
|
|
unsev(i)=exp(-ff(i)*ff(i)/eps**2) |
251 |
|
|
ENDDO |
252 |
|
|
ELSE |
253 |
|
|
DO i=1,ip1jm |
254 |
|
|
unsev(i)=eps**2/(ff(i)*ff(i)+eps**2) |
255 |
|
|
ENDDO |
256 |
|
|
END IF ! slab_sverdrup |
257 |
|
|
! 1/beta |
258 |
|
|
DO i=1,jjp1-1 |
259 |
|
|
unsbv((i-1)*iip1+1:i*iip1)=unsev((i-1)*iip1+1:i*iip1)/beta(i) |
260 |
|
|
END DO |
261 |
|
|
! 1/f |
262 |
|
|
ff=SIGN(MAX(ABS(ff),eps/100.),ff) ! avoid value 0 at equator... |
263 |
|
|
DO i=1,ip1jm |
264 |
|
|
unsfv(i)=(1.-unsev(i))/ff(i) |
265 |
|
|
END DO |
266 |
|
|
! compute values on 2D u grid |
267 |
|
|
! 1/eps |
268 |
|
|
unsev(:)=unsev(:)/eps |
269 |
|
|
CALL gr_v_scal(1,unsfv,unsfu) |
270 |
|
|
CALL gr_v_scal(1,unsev,unseu) |
271 |
|
|
END IF |
272 |
|
|
|
273 |
|
|
END SUBROUTINE ini_slab_transp |
274 |
|
|
|
275 |
|
|
SUBROUTINE divgrad_phy(nlevs,temp,delta) |
276 |
|
|
! Computes temperature tendency due to horizontal diffusion : |
277 |
|
|
! T Laplacian, later multiplied by diffusion coef and time-step |
278 |
|
|
|
279 |
|
|
IMPLICIT NONE |
280 |
|
|
|
281 |
|
|
INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers |
282 |
|
|
REAL, INTENT(IN) :: temp(klon_glo,nlevs) ! slab temperature |
283 |
|
|
REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.) |
284 |
|
|
REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs) |
285 |
|
|
REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs) |
286 |
|
|
INTEGER :: ll,iip1,jjp1 |
287 |
|
|
|
288 |
|
|
iip1=nbp_lon+1 |
289 |
|
|
jjp1=nbp_lat |
290 |
|
|
|
291 |
|
|
! transpose temp to 2D horiz. grid |
292 |
|
|
CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d) |
293 |
|
|
! computes gradient (proportional to heat flx) |
294 |
|
|
CALL grad(nlevs,delta_2d,ghx,ghy) |
295 |
|
|
! put flux to 0 at ocean / continent boundary |
296 |
|
|
DO ll=1,nlevs |
297 |
|
|
ghx(:,ll)=ghx(:,ll)*zmasqu |
298 |
|
|
ghy(:,ll)=ghy(:,ll)*zmasqv |
299 |
|
|
END DO |
300 |
|
|
! flux divergence |
301 |
|
|
CALL diverg(nlevs,ghx,ghy,delta_2d) |
302 |
|
|
! laplacian back to 1D grid |
303 |
|
|
CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta) |
304 |
|
|
|
305 |
|
|
RETURN |
306 |
|
|
END SUBROUTINE divgrad_phy |
307 |
|
|
|
308 |
|
|
SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy) |
309 |
|
|
! 1.5 Layer Ekman transport temperature tendency |
310 |
|
|
! note : tendency dt later multiplied by (delta t)/(rho.H) |
311 |
|
|
! to convert from divergence of heat fluxes to T |
312 |
|
|
|
313 |
|
|
IMPLICIT NONE |
314 |
|
|
|
315 |
|
|
! tx, ty : wind stress (different grids) |
316 |
|
|
! fluxm, fluz : mass *or heat* fluxes |
317 |
|
|
! dt : temperature tendency |
318 |
|
|
INTEGER ij |
319 |
|
|
|
320 |
|
|
! ts surface temp, td deep temp (diagnosed) |
321 |
|
|
REAL ts_phy(klon_glo) |
322 |
|
|
REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat) |
323 |
|
|
! zonal and meridional wind stress |
324 |
|
|
REAL tx_phy(klon_glo),ty_phy(klon_glo) |
325 |
|
|
REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat) |
326 |
|
|
REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1)) |
327 |
|
|
REAL tcurl((nbp_lon+1)*(nbp_lat-1)) |
328 |
|
|
! zonal and meridional Ekman mass fluxes at u, v points (2D grid) |
329 |
|
|
REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1)) |
330 |
|
|
! vertical and absolute mass fluxes (to estimate alpha) |
331 |
|
|
REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1)) |
332 |
|
|
! temperature tendency |
333 |
|
|
REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo) |
334 |
|
|
REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef |
335 |
|
|
|
336 |
|
|
INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
337 |
|
|
|
338 |
|
|
! Grid definitions |
339 |
|
|
iim=nbp_lon |
340 |
|
|
iip1=nbp_lon+1 |
341 |
|
|
iip2=nbp_lon+2 |
342 |
|
|
jjp1=nbp_lat |
343 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
344 |
|
|
ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
345 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
346 |
|
|
|
347 |
|
|
! Convert taux,y to 2D scalar grid |
348 |
|
|
! Note: 2D grid size = iim*jjm. iip1=iim+1 |
349 |
|
|
! First and last points in zonal direction are the same |
350 |
|
|
! we use 1 index ij from 1 to (iim+1)*(jjm+1) |
351 |
|
|
! north and south poles |
352 |
|
|
tx_phy(1)=0. |
353 |
|
|
tx_phy(klon_glo)=0. |
354 |
|
|
ty_phy(1)=0. |
355 |
|
|
ty_phy(klon_glo)=0. |
356 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu) |
357 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu) |
358 |
|
|
! convert to u,v grid (Arakawa C) |
359 |
|
|
! Multiply by f or eps to get mass flux |
360 |
|
|
! Meridional mass flux |
361 |
|
|
CALL gr_scal_v(1,txu,txv) ! wind stress at v points |
362 |
|
|
IF (slab_sverdrup) THEN ! Sverdrup bal. near equator |
363 |
|
|
tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:) |
364 |
|
|
fluxm=-tcurl*unsbv-txv*unsfv ! in kg.s-1.m-1 (zonal distance) |
365 |
|
|
ELSE |
366 |
|
|
CALL gr_scal_v(1,tyu,tyv) |
367 |
|
|
fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance) |
368 |
|
|
ENDIF |
369 |
|
|
! Zonal mass flux |
370 |
|
|
CALL gr_scal_u(1,txu,txu) ! wind stress at u points |
371 |
|
|
CALL gr_scal_u(1,tyu,tyu) |
372 |
|
|
fluxz=tyu*unsfu+txu*unseu |
373 |
|
|
|
374 |
|
|
! Correct flux: continent mask and horiz grid size |
375 |
|
|
! multiply m-flux by mask and dx: flux in kg.s-1 |
376 |
|
|
fluxm=fluxm*cv*cuvsurcv*zmasqv |
377 |
|
|
! multiply z-flux by mask and dy: flux in kg.s-1 |
378 |
|
|
fluxz=fluxz*cu*cvusurcu*zmasqu |
379 |
|
|
|
380 |
|
|
! Compute vertical and absolute mass flux (for variable alpha) |
381 |
|
|
IF (alpha_var) THEN |
382 |
|
|
DO ij=iip2,ip1jm |
383 |
|
|
fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
384 |
|
|
fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) & |
385 |
|
|
+ABS(fluxm(ij))+ABS(fluxm(ij-iip1)) |
386 |
|
|
ENDDO |
387 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
388 |
|
|
fluxt(ij+1)=fluxt(ij+iip1) |
389 |
|
|
fluxv(ij+1)=fluxv(ij+iip1) |
390 |
|
|
END DO |
391 |
|
|
fluxt(1)=SUM(ABS(fluxm(1:iim))) |
392 |
|
|
fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1))) |
393 |
|
|
fluxv(1)=-SUM(fluxm(1:iim)) |
394 |
|
|
fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
395 |
|
|
fluxt=MAX(fluxt,1.e-10) |
396 |
|
|
ENDIF |
397 |
|
|
|
398 |
|
|
! Compute alpha coefficient. |
399 |
|
|
! Tdeep = Tsurf * alpha + 271.35 * (1-alpha) |
400 |
|
|
IF (alpha_var) THEN |
401 |
|
|
! increase alpha (and Tdeep) in downwelling regions |
402 |
|
|
! and decrease in upwelling regions |
403 |
|
|
! to avoid "hot spots" where there is surface convergence |
404 |
|
|
DO ij=iip2,ip1jm |
405 |
|
|
alpha(ij)=alpham-fluxv(ij)/fluxt(ij)*(1.-alpham) |
406 |
|
|
ENDDO |
407 |
|
|
alpha(1:iip1)=alpham-fluxv(1)/fluxt(1)*(1.-alpham) |
408 |
|
|
alpha(ip1jm+1:ip1jmp1)=alpham-fluxv(ip1jmp1)/fluxt(ip1jmp1)*(1.-alpham) |
409 |
|
|
ELSE |
410 |
|
|
alpha(:)=alpham |
411 |
|
|
! Tsurf-Tdeep ~ 10� in the Tropics |
412 |
|
|
ENDIF |
413 |
|
|
|
414 |
|
|
! Estimate deep temperature |
415 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts) |
416 |
|
|
DO ij=1,ip1jmp1 |
417 |
|
|
td(ij)=271.35+(ts(ij)-271.35)*alpha(ij) |
418 |
|
|
td(ij)=MIN(td(ij),ts(ij)) |
419 |
|
|
END DO |
420 |
|
|
|
421 |
|
|
! Meridional heat flux: multiply mass flux by (ts-td) |
422 |
|
|
! flux in K.kg.s-1 |
423 |
|
|
IF (slab_upstream) THEN |
424 |
|
|
! upstream scheme to avoid hot spots |
425 |
|
|
DO ij=1,ip1jm |
426 |
|
|
IF (fluxm(ij).GE.0.) THEN |
427 |
|
|
fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij)) |
428 |
|
|
ELSE |
429 |
|
|
fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1)) |
430 |
|
|
END IF |
431 |
|
|
END DO |
432 |
|
|
ELSE |
433 |
|
|
! centered scheme better in mid-latitudes |
434 |
|
|
DO ij=1,ip1jm |
435 |
|
|
fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2. |
436 |
|
|
END DO |
437 |
|
|
ENDIF |
438 |
|
|
|
439 |
|
|
! Zonal heat flux |
440 |
|
|
! upstream scheme |
441 |
|
|
DO ij=iip2,ip1jm |
442 |
|
|
fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2. |
443 |
|
|
END DO |
444 |
|
|
DO ij=iip1*2,ip1jmp1,iip1 |
445 |
|
|
fluxz(ij)=fluxz(ij-iim) |
446 |
|
|
END DO |
447 |
|
|
|
448 |
|
|
! temperature tendency = divergence of heat fluxes |
449 |
|
|
! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface) |
450 |
|
|
DO ij=iip2,ip1jm |
451 |
|
|
dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) & |
452 |
|
|
/aire(ij) ! aire : grid area |
453 |
|
|
END DO |
454 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
455 |
|
|
dt(ij+1)=dt(ij+iip1) |
456 |
|
|
END DO |
457 |
|
|
! special treatment at the Poles |
458 |
|
|
dt(1)=SUM(fluxm(1:iim))/apoln |
459 |
|
|
dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols |
460 |
|
|
dt(2:iip1)=dt(1) |
461 |
|
|
dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1) |
462 |
|
|
|
463 |
|
|
! tendencies back to 1D grid |
464 |
|
|
CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy) |
465 |
|
|
|
466 |
|
|
RETURN |
467 |
|
|
END SUBROUTINE slab_ekman1 |
468 |
|
|
|
469 |
|
|
SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy_ek,dt_phy_gm,slab_gm) |
470 |
|
|
! Temperature tendency for 2-layers slab ocean |
471 |
|
|
! note : tendency dt later multiplied by (delta time)/(rho.H) |
472 |
|
|
! to convert from divergence of heat fluxes to T |
473 |
|
|
|
474 |
|
|
! 11/16 : Inclusion of GM-like eddy advection |
475 |
|
|
|
476 |
|
|
IMPLICIT NONE |
477 |
|
|
|
478 |
|
|
LOGICAL,INTENT(in) :: slab_gm |
479 |
|
|
! Here, temperature and flux variables are on 2 layers |
480 |
|
|
INTEGER ij |
481 |
|
|
|
482 |
|
|
! wind stress variables |
483 |
|
|
REAL tx_phy(klon_glo),ty_phy(klon_glo) |
484 |
|
|
REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1)) |
485 |
|
|
REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat) |
486 |
|
|
REAL tcurl((nbp_lon+1)*(nbp_lat-1)) |
487 |
|
|
! slab temperature on 1D, 2D grid |
488 |
|
|
REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2) |
489 |
|
|
! Temperature gradient, v-points |
490 |
|
|
REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat) |
491 |
|
|
! Vertical temperature difference, V-points |
492 |
|
|
REAL dtz((nbp_lon+1)*(nbp_lat-1)) |
493 |
|
|
! zonal and meridional mass fluxes at u, v points (2D grid) |
494 |
|
|
REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1)) |
495 |
|
|
! vertical mass flux between the 2 layers |
496 |
|
|
REAL fluxv_ek((nbp_lon+1)*nbp_lat) |
497 |
|
|
REAL fluxv_gm((nbp_lon+1)*nbp_lat) |
498 |
|
|
! zonal and meridional heat fluxes |
499 |
|
|
REAL fluxtz((nbp_lon+1)*nbp_lat,2) |
500 |
|
|
REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2) |
501 |
|
|
! temperature tendency (in K.s-1.kg.m-2) |
502 |
|
|
REAL dt_ek((nbp_lon+1)*nbp_lat,2), dt_phy_ek(klon_glo,2) |
503 |
|
|
REAL dt_gm((nbp_lon+1)*nbp_lat,2), dt_phy_gm(klon_glo,2) |
504 |
|
|
! helper vars |
505 |
|
|
REAL zonavg, fluxv |
506 |
|
|
REAL, PARAMETER :: sea_den=1025. ! sea water density |
507 |
|
|
|
508 |
|
|
INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
509 |
|
|
|
510 |
|
|
! Grid definitions |
511 |
|
|
iim=nbp_lon |
512 |
|
|
iip1=nbp_lon+1 |
513 |
|
|
iip2=nbp_lon+2 |
514 |
|
|
jjp1=nbp_lat |
515 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
516 |
|
|
ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
517 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
518 |
|
|
! Convert temperature to 2D grid |
519 |
|
|
CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts) |
520 |
|
|
|
521 |
|
|
! ------------------------------------ |
522 |
|
|
! Ekman mass fluxes and Temp tendency |
523 |
|
|
! ------------------------------------ |
524 |
|
|
! Convert taux,y to 2D scalar grid |
525 |
|
|
! north and south poles tx,ty no meaning |
526 |
|
|
tx_phy(1)=0. |
527 |
|
|
tx_phy(klon_glo)=0. |
528 |
|
|
ty_phy(1)=0. |
529 |
|
|
ty_phy(klon_glo)=0. |
530 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu) |
531 |
|
|
CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu) |
532 |
|
|
IF (ekman_zonavg) THEN ! use zonal average of wind stress |
533 |
|
|
DO ij=1,jjp1-2 |
534 |
|
|
zonavg=SUM(txu(ij*iip1+1:ij*iip1+iim))/iim |
535 |
|
|
txu(ij*iip1+1:(ij+1)*iip1)=zonavg |
536 |
|
|
zonavg=SUM(tyu(ij*iip1+1:ij*iip1+iim))/iim |
537 |
|
|
tyu(ij*iip1+1:(ij+1)*iip1)=zonavg |
538 |
|
|
END DO |
539 |
|
|
END IF |
540 |
|
|
|
541 |
|
|
! Divide taux,y by f or eps, and convert to 2D u,v grids |
542 |
|
|
! (Arakawa C grid) |
543 |
|
|
! Meridional flux |
544 |
|
|
CALL gr_scal_v(1,txu,txv) ! wind stress at v points |
545 |
|
|
fluxm=-txv*unsfv ! in kg.s-1.m-1 (zonal distance) |
546 |
|
|
IF (slab_sverdrup) THEN ! Sverdrup bal. near equator |
547 |
|
|
tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:) ! dtx/dy |
548 |
|
|
!poles curl = 0 |
549 |
|
|
tcurl(1:iip1)=0. |
550 |
|
|
tcurl(ip1jmi1+1:ip1jm)=0. |
551 |
|
|
fluxm=fluxm-tcurl*unsbv |
552 |
|
|
ENDIF |
553 |
|
|
IF (slab_tyeq) THEN ! meridional wind forcing at equator |
554 |
|
|
CALL gr_scal_v(1,tyu,tyv) |
555 |
|
|
fluxm=fluxm+tyv*unsev ! in kg.s-1.m-1 (zonal distance) |
556 |
|
|
ENDIF |
557 |
|
|
! apply continent mask, multiply by horiz grid dimension |
558 |
|
|
fluxm=fluxm*cv*cuvsurcv*zmasqv |
559 |
|
|
|
560 |
|
|
! Zonal flux |
561 |
|
|
IF (ekman_zonadv) THEN |
562 |
|
|
CALL gr_scal_u(1,txu,txu) ! wind stress at u points |
563 |
|
|
CALL gr_scal_u(1,tyu,tyu) |
564 |
|
|
fluxz=tyu*unsfu+txu*unseu |
565 |
|
|
! apply continent mask, multiply by horiz grid dimension |
566 |
|
|
fluxz=fluxz*cu*cvusurcu*zmasqu |
567 |
|
|
END IF |
568 |
|
|
|
569 |
|
|
! Vertical mass flux from mass budget (divergence of horiz fluxes) |
570 |
|
|
IF (ekman_zonadv) THEN |
571 |
|
|
DO ij=iip2,ip1jm |
572 |
|
|
fluxv_ek(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
573 |
|
|
ENDDO |
574 |
|
|
ELSE |
575 |
|
|
DO ij=iip2,ip1jm |
576 |
|
|
fluxv_ek(ij)=-fluxm(ij)+fluxm(ij-iip1) |
577 |
|
|
ENDDO |
578 |
|
|
END IF |
579 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
580 |
|
|
fluxv_ek(ij+1)=fluxv_ek(ij+iip1) |
581 |
|
|
END DO |
582 |
|
|
! vertical mass flux at Poles |
583 |
|
|
fluxv_ek(1)=-SUM(fluxm(1:iim)) |
584 |
|
|
fluxv_ek(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
585 |
|
|
|
586 |
|
|
! Meridional heat fluxes |
587 |
|
|
DO ij=1,ip1jm |
588 |
|
|
! centered scheme |
589 |
|
|
fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
590 |
|
|
fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
591 |
|
|
END DO |
592 |
|
|
|
593 |
|
|
! Zonal heat fluxes |
594 |
|
|
! Schema upstream |
595 |
|
|
IF (ekman_zonadv) THEN |
596 |
|
|
DO ij=iip2,ip1jm |
597 |
|
|
IF (fluxz(ij).GE.0.) THEN |
598 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
599 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
600 |
|
|
ELSE |
601 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
602 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
603 |
|
|
ENDIF |
604 |
|
|
ENDDO |
605 |
|
|
DO ij=iip1*2,ip1jmp1,iip1 |
606 |
|
|
fluxtz(ij,:)=fluxtz(ij-iim,:) |
607 |
|
|
END DO |
608 |
|
|
ELSE |
609 |
|
|
fluxtz(:,:)=0. |
610 |
|
|
ENDIF |
611 |
|
|
|
612 |
|
|
! Temperature tendency, horizontal advection: |
613 |
|
|
DO ij=iip2,ip1jm |
614 |
|
|
dt_ek(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
615 |
|
|
+fluxtm(ij,:)-fluxtm(ij-iip1,:) |
616 |
|
|
END DO |
617 |
|
|
! Poles |
618 |
|
|
dt_ek(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
619 |
|
|
dt_ek(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
620 |
|
|
|
621 |
|
|
! ------------------------------------ |
622 |
|
|
! GM mass fluxes and Temp tendency |
623 |
|
|
! ------------------------------------ |
624 |
|
|
IF (slab_gm) THEN |
625 |
|
|
! Vertical Temperature difference T1-T2 on v-grid points |
626 |
|
|
CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz) |
627 |
|
|
dtz(:)=MAX(dtz(:),0.25) |
628 |
|
|
! Horizontal Temperature differences |
629 |
|
|
CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty) |
630 |
|
|
! Meridional flux = -k.s (s=dyT/dzT) |
631 |
|
|
! Continent mask, multiply by dz/dy |
632 |
|
|
fluxm=dty/dtz*500.*cuvsurcv*zmasqv |
633 |
|
|
! slope limitation, multiply by kappa |
634 |
|
|
fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty) |
635 |
|
|
! convert to kg/s |
636 |
|
|
fluxm(:)=fluxm(:)*sea_den |
637 |
|
|
! Zonal flux = 0. (temporary) |
638 |
|
|
fluxz(:)=0. |
639 |
|
|
! Vertical mass flux from mass budget (divergence of horiz fluxes) |
640 |
|
|
DO ij=iip2,ip1jm |
641 |
|
|
fluxv_gm(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
642 |
|
|
ENDDO |
643 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
644 |
|
|
fluxv_gm(ij+1)=fluxv_gm(ij+iip1) |
645 |
|
|
END DO |
646 |
|
|
! vertical mass flux at Poles |
647 |
|
|
fluxv_gm(1)=-SUM(fluxm(1:iim)) |
648 |
|
|
fluxv_gm(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
649 |
|
|
|
650 |
|
|
! Meridional heat fluxes |
651 |
|
|
DO ij=1,ip1jm |
652 |
|
|
! centered scheme |
653 |
|
|
fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
654 |
|
|
fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
655 |
|
|
END DO |
656 |
|
|
|
657 |
|
|
! Zonal heat fluxes |
658 |
|
|
! Schema upstream |
659 |
|
|
DO ij=iip2,ip1jm |
660 |
|
|
IF (fluxz(ij).GE.0.) THEN |
661 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
662 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
663 |
|
|
ELSE |
664 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
665 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
666 |
|
|
ENDIF |
667 |
|
|
ENDDO |
668 |
|
|
DO ij=iip1*2,ip1jmp1,iip1 |
669 |
|
|
fluxtz(ij,:)=fluxtz(ij-iim,:) |
670 |
|
|
END DO |
671 |
|
|
|
672 |
|
|
! Temperature tendency : |
673 |
|
|
! divergence of horizontal heat fluxes |
674 |
|
|
DO ij=iip2,ip1jm |
675 |
|
|
dt_gm(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
676 |
|
|
+fluxtm(ij,:)-fluxtm(ij-iip1,:) |
677 |
|
|
END DO |
678 |
|
|
! Poles |
679 |
|
|
dt_gm(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
680 |
|
|
dt_gm(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
681 |
|
|
ELSE |
682 |
|
|
dt_gm(:,:)=0. |
683 |
|
|
fluxv_gm(:)=0. |
684 |
|
|
ENDIF ! slab_gm |
685 |
|
|
|
686 |
|
|
! ------------------------------------ |
687 |
|
|
! Temp tendency from vertical advection |
688 |
|
|
! Divide by cell area |
689 |
|
|
! ------------------------------------ |
690 |
|
|
! vertical heat flux = mass flux * T, upstream scheme |
691 |
|
|
DO ij=iip2,ip1jm |
692 |
|
|
fluxv=fluxv_ek(ij)+fluxv_gm(ij) ! net flux, needed for upstream scheme |
693 |
|
|
IF (fluxv.GT.0.) THEN |
694 |
|
|
dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,2) |
695 |
|
|
dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,2) |
696 |
|
|
dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,2) |
697 |
|
|
dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,2) |
698 |
|
|
ELSE |
699 |
|
|
dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,1) |
700 |
|
|
dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,1) |
701 |
|
|
dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,1) |
702 |
|
|
dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,1) |
703 |
|
|
ENDIF |
704 |
|
|
! divide by cell area |
705 |
|
|
dt_ek(ij,:)=dt_ek(ij,:)/aire(ij) |
706 |
|
|
dt_gm(ij,:)=dt_gm(ij,:)/aire(ij) |
707 |
|
|
END DO |
708 |
|
|
! North Pole |
709 |
|
|
fluxv=fluxv_ek(1)+fluxv_gm(1) |
710 |
|
|
IF (fluxv.GT.0.) THEN |
711 |
|
|
dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,2) |
712 |
|
|
dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,2) |
713 |
|
|
dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,2) |
714 |
|
|
dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,2) |
715 |
|
|
ELSE |
716 |
|
|
dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,1) |
717 |
|
|
dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,1) |
718 |
|
|
dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,1) |
719 |
|
|
dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,1) |
720 |
|
|
ENDIF |
721 |
|
|
dt_ek(1,:)=dt_ek(1,:)/apoln |
722 |
|
|
dt_gm(1,:)=dt_gm(1,:)/apoln |
723 |
|
|
! South pole |
724 |
|
|
fluxv=fluxv_ek(ip1jmp1)+fluxv_gm(ip1jmp1) |
725 |
|
|
IF (fluxv.GT.0.) THEN |
726 |
|
|
dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,2) |
727 |
|
|
dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,2) |
728 |
|
|
dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,2) |
729 |
|
|
dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,2) |
730 |
|
|
ELSE |
731 |
|
|
dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,1) |
732 |
|
|
dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,1) |
733 |
|
|
dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,1) |
734 |
|
|
dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,1) |
735 |
|
|
ENDIF |
736 |
|
|
dt_ek(ip1jmp1,:)=dt_ek(ip1jmp1,:)/apols |
737 |
|
|
dt_gm(ip1jmp1,:)=dt_gm(ip1jmp1,:)/apols |
738 |
|
|
|
739 |
|
|
dt_ek(2:iip1,1)=dt_ek(1,1) |
740 |
|
|
dt_ek(2:iip1,2)=dt_ek(1,2) |
741 |
|
|
dt_gm(2:iip1,1)=dt_gm(1,1) |
742 |
|
|
dt_gm(2:iip1,2)=dt_gm(1,2) |
743 |
|
|
dt_ek(ip1jm+1:ip1jmp1,1)=dt_ek(ip1jmp1,1) |
744 |
|
|
dt_ek(ip1jm+1:ip1jmp1,2)=dt_ek(ip1jmp1,2) |
745 |
|
|
dt_gm(ip1jm+1:ip1jmp1,1)=dt_gm(ip1jmp1,1) |
746 |
|
|
dt_gm(ip1jm+1:ip1jmp1,2)=dt_gm(ip1jmp1,2) |
747 |
|
|
|
748 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
749 |
|
|
dt_gm(ij+1,:)=dt_gm(ij+iip1,:) |
750 |
|
|
dt_ek(ij+1,:)=dt_ek(ij+iip1,:) |
751 |
|
|
END DO |
752 |
|
|
|
753 |
|
|
! T tendency back to 1D grid... |
754 |
|
|
CALL gr_dyn_fi(2,iip1,jjp1,dt_ek,dt_phy_ek) |
755 |
|
|
CALL gr_dyn_fi(2,iip1,jjp1,dt_gm,dt_phy_gm) |
756 |
|
|
|
757 |
|
|
RETURN |
758 |
|
|
END SUBROUTINE slab_ekman2 |
759 |
|
|
|
760 |
|
|
SUBROUTINE slab_gmdiff(ts_phy,dt_phy) |
761 |
|
|
! Temperature tendency for 2-layers slab ocean |
762 |
|
|
! Due to Gent-McWilliams type eddy-induced advection |
763 |
|
|
|
764 |
|
|
IMPLICIT NONE |
765 |
|
|
|
766 |
|
|
! Here, temperature and flux variables are on 2 layers |
767 |
|
|
INTEGER ij |
768 |
|
|
! Temperature gradient, v-points |
769 |
|
|
REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat) |
770 |
|
|
! Vertical temperature difference, V-points |
771 |
|
|
REAL dtz((nbp_lon+1)*(nbp_lat-1)) |
772 |
|
|
! slab temperature on 1D, 2D grid |
773 |
|
|
REAL ts_phy(klon_glo,2),ts((nbp_lon+1)*nbp_lat,2) |
774 |
|
|
! zonal and meridional mass fluxes at u, v points (2D grid) |
775 |
|
|
REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1)) |
776 |
|
|
! vertical mass flux between the 2 layers |
777 |
|
|
REAL fluxv((nbp_lon+1)*nbp_lat) |
778 |
|
|
! zonal and meridional heat fluxes |
779 |
|
|
REAL fluxtz((nbp_lon+1)*nbp_lat,2) |
780 |
|
|
REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2) |
781 |
|
|
! temperature tendency (in K.s-1.kg.m-2) |
782 |
|
|
REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2) |
783 |
|
|
|
784 |
|
|
INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
785 |
|
|
|
786 |
|
|
! Grid definitions |
787 |
|
|
iim=nbp_lon |
788 |
|
|
iip1=nbp_lon+1 |
789 |
|
|
iip2=nbp_lon+2 |
790 |
|
|
jjp1=nbp_lat |
791 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
792 |
|
|
ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
793 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
794 |
|
|
|
795 |
|
|
! Convert temperature to 2D grid |
796 |
|
|
CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts) |
797 |
|
|
! Vertical Temperature difference T1-T2 on v-grid points |
798 |
|
|
CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz) |
799 |
|
|
dtz(:)=MAX(dtz(:),0.25) |
800 |
|
|
! Horizontal Temperature differences |
801 |
|
|
CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty) |
802 |
|
|
! Meridional flux = -k.s (s=dyT/dzT) |
803 |
|
|
! Continent mask, multiply by dz/dy |
804 |
|
|
fluxm=dty/dtz*500.*cuvsurcv*zmasqv |
805 |
|
|
! slope limitation, multiply by kappa |
806 |
|
|
fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty) |
807 |
|
|
! Zonal flux = 0. (temporary) |
808 |
|
|
fluxz(:)=0. |
809 |
|
|
! Vertical mass flux from mass budget (divergence of horiz fluxes) |
810 |
|
|
DO ij=iip2,ip1jm |
811 |
|
|
fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
812 |
|
|
ENDDO |
813 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
814 |
|
|
fluxv(ij+1)=fluxv(ij+iip1) |
815 |
|
|
END DO |
816 |
|
|
! vertical mass flux at Poles |
817 |
|
|
fluxv(1)=-SUM(fluxm(1:iim)) |
818 |
|
|
fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
819 |
|
|
fluxv=fluxv |
820 |
|
|
|
821 |
|
|
! Meridional heat fluxes |
822 |
|
|
DO ij=1,ip1jm |
823 |
|
|
! centered scheme |
824 |
|
|
fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
825 |
|
|
fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
826 |
|
|
END DO |
827 |
|
|
|
828 |
|
|
! Zonal heat fluxes |
829 |
|
|
! Schema upstream |
830 |
|
|
DO ij=iip2,ip1jm |
831 |
|
|
IF (fluxz(ij).GE.0.) THEN |
832 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
833 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
834 |
|
|
ELSE |
835 |
|
|
fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
836 |
|
|
fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
837 |
|
|
ENDIF |
838 |
|
|
ENDDO |
839 |
|
|
DO ij=iip1*2,ip1jmp1,iip1 |
840 |
|
|
fluxtz(ij,:)=fluxtz(ij-iim,:) |
841 |
|
|
END DO |
842 |
|
|
|
843 |
|
|
! Temperature tendency : |
844 |
|
|
DO ij=iip2,ip1jm |
845 |
|
|
! divergence of horizontal heat fluxes |
846 |
|
|
dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
847 |
|
|
+fluxtm(ij,:)-fluxtm(ij-iip1,:) |
848 |
|
|
! + vertical heat flux (mass flux * T, upstream scheme) |
849 |
|
|
IF (fluxv(ij).GT.0.) THEN |
850 |
|
|
dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2) |
851 |
|
|
dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2) |
852 |
|
|
ELSE |
853 |
|
|
dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1) |
854 |
|
|
dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1) |
855 |
|
|
ENDIF |
856 |
|
|
! divide by cell area |
857 |
|
|
dt(ij,:)=dt(ij,:)/aire(ij) |
858 |
|
|
END DO |
859 |
|
|
DO ij=iip1,ip1jmi1,iip1 |
860 |
|
|
dt(ij+1,:)=dt(ij+iip1,:) |
861 |
|
|
END DO |
862 |
|
|
! Poles |
863 |
|
|
dt(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
864 |
|
|
IF (fluxv(1).GT.0.) THEN |
865 |
|
|
dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2) |
866 |
|
|
dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2) |
867 |
|
|
ELSE |
868 |
|
|
dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1) |
869 |
|
|
dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1) |
870 |
|
|
ENDIF |
871 |
|
|
dt(1,:)=dt(1,:)/apoln |
872 |
|
|
dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
873 |
|
|
IF (fluxv(ip1jmp1).GT.0.) THEN |
874 |
|
|
dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2) |
875 |
|
|
dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2) |
876 |
|
|
ELSE |
877 |
|
|
dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1) |
878 |
|
|
dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1) |
879 |
|
|
ENDIF |
880 |
|
|
dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols |
881 |
|
|
dt(2:iip1,1)=dt(1,1) |
882 |
|
|
dt(2:iip1,2)=dt(1,2) |
883 |
|
|
dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1) |
884 |
|
|
dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2) |
885 |
|
|
|
886 |
|
|
! T tendency back to 1D grid... |
887 |
|
|
CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy) |
888 |
|
|
|
889 |
|
|
RETURN |
890 |
|
|
END SUBROUTINE slab_gmdiff |
891 |
|
|
|
892 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
893 |
|
|
|
894 |
|
|
SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn) |
895 |
|
|
! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid |
896 |
|
|
IMPLICIT NONE |
897 |
|
|
|
898 |
|
|
INTEGER,INTENT(IN) :: im,jm,nfield |
899 |
|
|
REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid |
900 |
|
|
REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid |
901 |
|
|
|
902 |
|
|
INTEGER :: i,j,ifield,ig |
903 |
|
|
|
904 |
|
|
DO ifield=1,nfield |
905 |
|
|
! Handle poles |
906 |
|
|
DO i=1,im |
907 |
|
|
pdyn(i,1,ifield)=pfi(1,ifield) |
908 |
|
|
pdyn(i,jm,ifield)=pfi(klon_glo,ifield) |
909 |
|
|
ENDDO |
910 |
|
|
! Other points |
911 |
|
|
DO j=2,jm-1 |
912 |
|
|
ig=2+(j-2)*(im-1) |
913 |
|
|
CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1) |
914 |
|
|
pdyn(im,j,ifield)=pdyn(1,j,ifield) |
915 |
|
|
ENDDO |
916 |
|
|
ENDDO ! of DO ifield=1,nfield |
917 |
|
|
|
918 |
|
|
END SUBROUTINE gr_fi_dyn |
919 |
|
|
|
920 |
|
|
SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi) |
921 |
|
|
! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid |
922 |
|
|
IMPLICIT NONE |
923 |
|
|
|
924 |
|
|
INTEGER,INTENT(IN) :: im,jm,nfield |
925 |
|
|
REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid |
926 |
|
|
REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid |
927 |
|
|
|
928 |
|
|
INTEGER j,ifield,ig |
929 |
|
|
|
930 |
|
|
CHARACTER (len = 20) :: modname = 'slab_heat_transp' |
931 |
|
|
CHARACTER (len = 80) :: abort_message |
932 |
|
|
|
933 |
|
|
! Sanity check: |
934 |
|
|
IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN |
935 |
|
|
abort_message="gr_dyn_fi error, wrong sizes" |
936 |
|
|
CALL abort_physic(modname,abort_message,1) |
937 |
|
|
ENDIF |
938 |
|
|
|
939 |
|
|
! Handle poles |
940 |
|
|
CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo) |
941 |
|
|
CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo) |
942 |
|
|
! Other points |
943 |
|
|
DO ifield=1,nfield |
944 |
|
|
DO j=2,jm-1 |
945 |
|
|
ig=2+(j-2)*(im-1) |
946 |
|
|
CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1) |
947 |
|
|
ENDDO |
948 |
|
|
ENDDO |
949 |
|
|
|
950 |
|
|
END SUBROUTINE gr_dyn_fi |
951 |
|
|
|
952 |
|
|
SUBROUTINE grad(klevel,pg,pgx,pgy) |
953 |
|
|
! compute the covariant components pgx,pgy of the gradient of pg |
954 |
|
|
! pgx = d(pg)/dx * delta(x) = delta(pg) |
955 |
|
|
IMPLICIT NONE |
956 |
|
|
|
957 |
|
|
INTEGER,INTENT(IN) :: klevel |
958 |
|
|
REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel) |
959 |
|
|
REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel) |
960 |
|
|
REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel) |
961 |
|
|
|
962 |
|
|
INTEGER :: l,ij |
963 |
|
|
INTEGER :: iim,iip1,ip1jm,ip1jmp1 |
964 |
|
|
|
965 |
|
|
iim=nbp_lon |
966 |
|
|
iip1=nbp_lon+1 |
967 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
968 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
969 |
|
|
|
970 |
|
|
DO l=1,klevel |
971 |
|
|
DO ij=1,ip1jmp1-1 |
972 |
|
|
pgx(ij,l)=pg(ij+1,l)-pg(ij,l) |
973 |
|
|
ENDDO |
974 |
|
|
! correction for pgx(ip1,j,l) ... |
975 |
|
|
! ... pgx(iip1,j,l)=pgx(1,j,l) ... |
976 |
|
|
DO ij=iip1,ip1jmp1,iip1 |
977 |
|
|
pgx(ij,l)=pgx(ij-iim,l) |
978 |
|
|
ENDDO |
979 |
|
|
DO ij=1,ip1jm |
980 |
|
|
pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l) |
981 |
|
|
ENDDO |
982 |
|
|
ENDDO |
983 |
|
|
|
984 |
|
|
END SUBROUTINE grad |
985 |
|
|
|
986 |
|
|
SUBROUTINE diverg(klevel,x,y,div) |
987 |
|
|
! computes the divergence of a vector field of components |
988 |
|
|
! x,y. x and y being covariant components |
989 |
|
|
IMPLICIT NONE |
990 |
|
|
|
991 |
|
|
INTEGER,INTENT(IN) :: klevel |
992 |
|
|
REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel) |
993 |
|
|
REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel) |
994 |
|
|
REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel) |
995 |
|
|
|
996 |
|
|
INTEGER :: l,ij |
997 |
|
|
INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1 |
998 |
|
|
|
999 |
|
|
REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1) |
1000 |
|
|
REAL :: sumypn,sumyps |
1001 |
|
|
REAL,EXTERNAL :: SSUM |
1002 |
|
|
|
1003 |
|
|
iim=nbp_lon |
1004 |
|
|
iip1=nbp_lon+1 |
1005 |
|
|
iip2=nbp_lon+2 |
1006 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
1007 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
1008 |
|
|
ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
1009 |
|
|
|
1010 |
|
|
DO l=1,klevel |
1011 |
|
|
DO ij=iip2,ip1jm-1 |
1012 |
|
|
div(ij+1,l)= & |
1013 |
|
|
cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ & |
1014 |
|
|
cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l) |
1015 |
|
|
ENDDO |
1016 |
|
|
! correction for div(1,j,l) ... |
1017 |
|
|
! ... div(1,j,l)= div(iip1,j,l) ... |
1018 |
|
|
DO ij=iip2,ip1jm,iip1 |
1019 |
|
|
div(ij,l)=div(ij+iim,l) |
1020 |
|
|
ENDDO |
1021 |
|
|
! at the poles |
1022 |
|
|
DO ij=1,iim |
1023 |
|
|
aiy1(ij)=cuvsurcv(ij)*y(ij,l) |
1024 |
|
|
aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l) |
1025 |
|
|
ENDDO |
1026 |
|
|
sumypn=SSUM(iim,aiy1,1)/apoln |
1027 |
|
|
sumyps=SSUM(iim,aiy2,1)/apols |
1028 |
|
|
DO ij=1,iip1 |
1029 |
|
|
div(ij,l)=-sumypn |
1030 |
|
|
div(ij+ip1jm,l)=sumyps |
1031 |
|
|
ENDDO |
1032 |
|
|
! End (poles) |
1033 |
|
|
ENDDO ! of DO l=1,klevel |
1034 |
|
|
|
1035 |
|
|
!!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 ) |
1036 |
|
|
DO l=1,klevel |
1037 |
|
|
DO ij=iip2,ip1jm |
1038 |
|
|
div(ij,l)=div(ij,l)*unsaire(ij) |
1039 |
|
|
ENDDO |
1040 |
|
|
ENDDO |
1041 |
|
|
|
1042 |
|
|
END SUBROUTINE diverg |
1043 |
|
|
|
1044 |
|
|
SUBROUTINE gr_v_scal(nx,x_v,x_scal) |
1045 |
|
|
! convert values from v points to scalar points on C-grid |
1046 |
|
|
! used to compute unsfu, unseu (u points, but depends only on latitude) |
1047 |
|
|
IMPLICIT NONE |
1048 |
|
|
|
1049 |
|
|
INTEGER,INTENT(IN) :: nx ! number of levels or fields |
1050 |
|
|
REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
1051 |
|
|
REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
1052 |
|
|
|
1053 |
|
|
INTEGER :: l,ij |
1054 |
|
|
INTEGER :: iip1,iip2,ip1jm,ip1jmp1 |
1055 |
|
|
|
1056 |
|
|
iip1=nbp_lon+1 |
1057 |
|
|
iip2=nbp_lon+2 |
1058 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
1059 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
1060 |
|
|
|
1061 |
|
|
DO l=1,nx |
1062 |
|
|
DO ij=iip2,ip1jm |
1063 |
|
|
x_scal(ij,l)= & |
1064 |
|
|
(airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) & |
1065 |
|
|
/(airev(ij-iip1)+airev(ij)) |
1066 |
|
|
ENDDO |
1067 |
|
|
DO ij=1,iip1 |
1068 |
|
|
x_scal(ij,l)=0. |
1069 |
|
|
ENDDO |
1070 |
|
|
DO ij=ip1jm+1,ip1jmp1 |
1071 |
|
|
x_scal(ij,l)=0. |
1072 |
|
|
ENDDO |
1073 |
|
|
ENDDO |
1074 |
|
|
|
1075 |
|
|
END SUBROUTINE gr_v_scal |
1076 |
|
|
|
1077 |
|
|
SUBROUTINE gr_scal_v(nx,x_scal,x_v) |
1078 |
|
|
! convert values from scalar points to v points on C-grid |
1079 |
|
|
! used to compute wind stress at V points |
1080 |
|
|
IMPLICIT NONE |
1081 |
|
|
|
1082 |
|
|
INTEGER,INTENT(IN) :: nx ! number of levels or fields |
1083 |
|
|
REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
1084 |
|
|
REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
1085 |
|
|
|
1086 |
|
|
INTEGER :: l,ij |
1087 |
|
|
INTEGER :: iip1,ip1jm |
1088 |
|
|
|
1089 |
|
|
iip1=nbp_lon+1 |
1090 |
|
|
ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
1091 |
|
|
|
1092 |
|
|
DO l=1,nx |
1093 |
|
|
DO ij=1,ip1jm |
1094 |
|
|
x_v(ij,l)= & |
1095 |
|
|
(cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ & |
1096 |
|
|
cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) & |
1097 |
|
|
/(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1)) |
1098 |
|
|
ENDDO |
1099 |
|
|
ENDDO |
1100 |
|
|
|
1101 |
|
|
END SUBROUTINE gr_scal_v |
1102 |
|
|
|
1103 |
|
|
SUBROUTINE gr_scal_u(nx,x_scal,x_u) |
1104 |
|
|
! convert values from scalar points to U points on C-grid |
1105 |
|
|
! used to compute wind stress at U points |
1106 |
|
|
IMPLICIT NONE |
1107 |
|
|
|
1108 |
|
|
INTEGER,INTENT(IN) :: nx |
1109 |
|
|
REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx) |
1110 |
|
|
REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
1111 |
|
|
|
1112 |
|
|
INTEGER :: l,ij |
1113 |
|
|
INTEGER :: iip1,jjp1,ip1jmp1 |
1114 |
|
|
|
1115 |
|
|
iip1=nbp_lon+1 |
1116 |
|
|
jjp1=nbp_lat |
1117 |
|
|
ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
1118 |
|
|
|
1119 |
|
|
DO l=1,nx |
1120 |
|
|
DO ij=1,ip1jmp1-1 |
1121 |
|
|
x_u(ij,l)= & |
1122 |
|
|
(aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) & |
1123 |
|
|
/(aire(ij)+aire(ij+1)) |
1124 |
|
|
ENDDO |
1125 |
|
|
ENDDO |
1126 |
|
|
|
1127 |
|
|
CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1) |
1128 |
|
|
|
1129 |
|
|
END SUBROUTINE gr_scal_u |
1130 |
|
|
|
1131 |
|
|
END MODULE slab_heat_transp_mod |