source: LMDZ6/branches/Amaury_dev/libf/phylmd/slab_heat_transp_mod.F90

Last change on this file was 5160, checked in by abarral, 3 months ago

Put .h into modules
File size: 38.7 KB
Line 
1MODULE slab_heat_transp_mod
2
3  ! Slab ocean : temperature tendencies due to horizontal diffusion
4  ! and / or Ekman transport
5
6  USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat, klon_glo
7  USE lmdz_abort_physic, ONLY: abort_physic
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
80CONTAINS
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/=((nbp_lon + 1) * (nbp_lat - 1))).OR. &
113            (ip1jmp1/=((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 lmdz_ioipsl_getin_p, 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)>1e-5 .OR. zmasq_2d(i + 1)>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)>1e-5 .OR. zmasq_2d(i + iip1)>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>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  END SUBROUTINE divgrad_phy
306
307  SUBROUTINE slab_ekman1(tx_phy, ty_phy, ts_phy, dt_phy)
308    ! 1.5 Layer Ekman transport temperature tendency
309    ! note : tendency dt later multiplied by (delta t)/(rho.H)
310    ! to convert from divergence of heat fluxes to T
311
312    IMPLICIT NONE
313
314    ! tx, ty : wind stress (different grids)
315    ! fluxm, fluz : mass *or heat* fluxes
316    ! dt : temperature tendency
317    INTEGER ij
318
319    ! ts surface temp, td deep temp (diagnosed)
320    REAL ts_phy(klon_glo)
321    REAL ts((nbp_lon + 1) * nbp_lat), td((nbp_lon + 1) * nbp_lat)
322    ! zonal and meridional wind stress
323    REAL tx_phy(klon_glo), ty_phy(klon_glo)
324    REAL tyu((nbp_lon + 1) * nbp_lat), txu((nbp_lon + 1) * nbp_lat)
325    REAL txv((nbp_lon + 1) * (nbp_lat - 1)), tyv((nbp_lon + 1) * (nbp_lat - 1))
326    REAL tcurl((nbp_lon + 1) * (nbp_lat - 1))
327    ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
328    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
329    ! vertical  and absolute mass fluxes (to estimate alpha)
330    REAL fluxv((nbp_lon + 1) * nbp_lat), fluxt((nbp_lon + 1) * (nbp_lat - 1))
331    ! temperature tendency
332    REAL dt((nbp_lon + 1) * nbp_lat), dt_phy(klon_glo)
333    REAL alpha((nbp_lon + 1) * nbp_lat) ! deep temperature coef
334
335    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1
336
337    ! Grid definitions
338    iim = nbp_lon
339    iip1 = nbp_lon + 1
340    iip2 = nbp_lon + 2
341    jjp1 = nbp_lat
342    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
343    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
344    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
345
346    ! Convert taux,y to 2D  scalar grid
347    ! Note: 2D grid size = iim*jjm. iip1=iim+1
348    ! First and last points in zonal direction are the same
349    ! we use 1 index ij from 1 to (iim+1)*(jjm+1)
350    ! north and south poles
351    tx_phy(1) = 0.
352    tx_phy(klon_glo) = 0.
353    ty_phy(1) = 0.
354    ty_phy(klon_glo) = 0.
355    CALL gr_fi_dyn(1, iip1, jjp1, tx_phy, txu)
356    CALL gr_fi_dyn(1, iip1, jjp1, ty_phy, tyu)
357    ! convert to u,v grid (Arakawa C)
358    ! Multiply by f or eps to get mass flux
359    ! Meridional mass flux
360    CALL gr_scal_v(1, txu, txv) ! wind stress at v points
361    IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
362      tcurl = (txu(1:ip1jm) - txu(iip2:ip1jmp1)) / cv(:)
363      fluxm = -tcurl * unsbv - txv * unsfv ! in kg.s-1.m-1 (zonal distance)
364    ELSE
365      CALL gr_scal_v(1, tyu, tyv)
366      fluxm = tyv * unsev - txv * unsfv ! in kg.s-1.m-1 (zonal distance)
367    ENDIF
368    ! Zonal mass flux
369    CALL gr_scal_u(1, txu, txu) ! wind stress at u points
370    CALL gr_scal_u(1, tyu, tyu)
371    fluxz = tyu * unsfu + txu * unseu
372
373    ! Correct flux: continent mask and horiz grid size
374    ! multiply m-flux by mask and dx: flux in kg.s-1
375    fluxm = fluxm * cv * cuvsurcv * zmasqv
376    ! multiply z-flux by mask and dy: flux in kg.s-1
377    fluxz = fluxz * cu * cvusurcu * zmasqu
378
379    ! Compute vertical  and absolute mass flux (for variable alpha)
380    IF (alpha_var) THEN
381      DO ij = iip2, ip1jm
382        fluxv(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
383        fluxt(ij) = ABS(fluxz(ij)) + ABS(fluxz(ij - 1)) &
384                + ABS(fluxm(ij)) + ABS(fluxm(ij - iip1))
385      ENDDO
386      DO ij = iip1, ip1jmi1, iip1
387        fluxt(ij + 1) = fluxt(ij + iip1)
388        fluxv(ij + 1) = fluxv(ij + iip1)
389      END DO
390      fluxt(1) = SUM(ABS(fluxm(1:iim)))
391      fluxt(ip1jmp1) = SUM(ABS(fluxm(ip1jm - iim:ip1jm - 1)))
392      fluxv(1) = -SUM(fluxm(1:iim))
393      fluxv(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
394      fluxt = MAX(fluxt, 1.e-10)
395    ENDIF
396
397    ! Compute alpha coefficient.
398    ! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
399    IF (alpha_var) THEN
400      ! increase alpha (and Tdeep) in downwelling regions
401      ! and decrease in upwelling regions
402      ! to avoid "hot spots" where there is surface convergence
403      DO ij = iip2, ip1jm
404        alpha(ij) = alpham - fluxv(ij) / fluxt(ij) * (1. - alpham)
405      ENDDO
406      alpha(1:iip1) = alpham - fluxv(1) / fluxt(1) * (1. - alpham)
407      alpha(ip1jm + 1:ip1jmp1) = alpham - fluxv(ip1jmp1) / fluxt(ip1jmp1) * (1. - alpham)
408    ELSE
409      alpha(:) = alpham
410      ! Tsurf-Tdeep ~ 10° in the Tropics
411    ENDIF
412
413    ! Estimate deep temperature
414    CALL gr_fi_dyn(1, iip1, jjp1, ts_phy, ts)
415    DO ij = 1, ip1jmp1
416      td(ij) = 271.35 + (ts(ij) - 271.35) * alpha(ij)
417      td(ij) = MIN(td(ij), ts(ij))
418    END DO
419
420    ! Meridional heat flux: multiply mass flux by (ts-td)
421    ! flux in K.kg.s-1
422    IF (slab_upstream) THEN
423      ! upstream scheme to avoid hot spots
424      DO ij = 1, ip1jm
425        IF (fluxm(ij)>=0.) THEN
426          fluxm(ij) = fluxm(ij) * (ts(ij + iip1) - td(ij))
427        ELSE
428          fluxm(ij) = fluxm(ij) * (ts(ij) - td(ij + iip1))
429        END IF
430      END DO
431    ELSE
432      ! centered scheme better in mid-latitudes
433      DO ij = 1, ip1jm
434        fluxm(ij) = fluxm(ij) * (ts(ij + iip1) + ts(ij) - td(ij) - td(ij + iip1)) / 2.
435      END DO
436    ENDIF
437
438    ! Zonal heat flux
439    ! upstream scheme
440    DO ij = iip2, ip1jm
441      fluxz(ij) = fluxz(ij) * (ts(ij) + ts(ij + 1) - td(ij + 1) - td(ij)) / 2.
442    END DO
443    DO ij = iip1 * 2, ip1jmp1, iip1
444      fluxz(ij) = fluxz(ij - iim)
445    END DO
446
447    ! temperature tendency = divergence of heat fluxes
448    ! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
449    DO ij = iip2, ip1jm
450      dt(ij) = (fluxz(ij - 1) - fluxz(ij) + fluxm(ij) - fluxm(ij - iip1)) &
451              / aire(ij) ! aire : grid area
452    END DO
453    DO ij = iip1, ip1jmi1, iip1
454      dt(ij + 1) = dt(ij + iip1)
455    END DO
456    ! special treatment at the Poles
457    dt(1) = SUM(fluxm(1:iim)) / apoln
458    dt(ip1jmp1) = -SUM(fluxm(ip1jm - iim:ip1jm - 1)) / apols
459    dt(2:iip1) = dt(1)
460    dt(ip1jm + 1:ip1jmp1) = dt(ip1jmp1)
461
462    ! tendencies back to 1D grid
463    CALL gr_dyn_fi(1, iip1, jjp1, dt, dt_phy)
464
465  END SUBROUTINE slab_ekman1
466
467  SUBROUTINE slab_ekman2(tx_phy, ty_phy, ts_phy, dt_phy_ek, dt_phy_gm, slab_gm)
468    ! Temperature tendency for 2-layers slab ocean
469    ! note : tendency dt later multiplied by (delta time)/(rho.H)
470    ! to convert from divergence of heat fluxes to T
471
472    ! 11/16 : Inclusion of GM-like eddy advection
473
474    IMPLICIT NONE
475
476    LOGICAL, INTENT(IN) :: slab_gm
477    ! Here, temperature and flux variables are on 2 layers
478    INTEGER ij
479
480    ! wind stress variables
481    REAL tx_phy(klon_glo), ty_phy(klon_glo)
482    REAL txv((nbp_lon + 1) * (nbp_lat - 1)), tyv((nbp_lon + 1) * (nbp_lat - 1))
483    REAL tyu((nbp_lon + 1) * nbp_lat), txu((nbp_lon + 1) * nbp_lat)
484    REAL tcurl((nbp_lon + 1) * (nbp_lat - 1))
485    ! slab temperature on  1D, 2D grid
486    REAL ts_phy(klon_glo, 2), ts((nbp_lon + 1) * nbp_lat, 2)
487    ! Temperature gradient, v-points
488    REAL dty((nbp_lon + 1) * (nbp_lat - 1)), dtx((nbp_lon + 1) * nbp_lat)
489    ! Vertical temperature difference, V-points
490    REAL dtz((nbp_lon + 1) * (nbp_lat - 1))
491    ! zonal and meridional mass fluxes at u, v points (2D grid)
492    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
493    ! vertical mass flux between the 2 layers
494    REAL fluxv_ek((nbp_lon + 1) * nbp_lat)
495    REAL fluxv_gm((nbp_lon + 1) * nbp_lat)
496    ! zonal and meridional heat fluxes
497    REAL fluxtz((nbp_lon + 1) * nbp_lat, 2)
498    REAL fluxtm((nbp_lon + 1) * (nbp_lat - 1), 2)
499    ! temperature tendency (in K.s-1.kg.m-2)
500    REAL dt_ek((nbp_lon + 1) * nbp_lat, 2), dt_phy_ek(klon_glo, 2)
501    REAL dt_gm((nbp_lon + 1) * nbp_lat, 2), dt_phy_gm(klon_glo, 2)
502    ! helper vars
503    REAL zonavg, fluxv
504    REAL, PARAMETER :: sea_den = 1025. ! sea water density
505
506    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1
507
508    ! Grid definitions
509    iim = nbp_lon
510    iip1 = nbp_lon + 1
511    iip2 = nbp_lon + 2
512    jjp1 = nbp_lat
513    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
514    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
515    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
516    ! Convert temperature to 2D grid
517    CALL gr_fi_dyn(2, iip1, jjp1, ts_phy, ts)
518
519    ! ------------------------------------
520    ! Ekman mass fluxes and Temp tendency
521    ! ------------------------------------
522    ! Convert taux,y to 2D  scalar grid
523    ! north and south poles tx,ty no meaning
524    tx_phy(1) = 0.
525    tx_phy(klon_glo) = 0.
526    ty_phy(1) = 0.
527    ty_phy(klon_glo) = 0.
528    CALL gr_fi_dyn(1, iip1, jjp1, tx_phy, txu)
529    CALL gr_fi_dyn(1, iip1, jjp1, ty_phy, tyu)
530    IF (ekman_zonavg) THEN ! use zonal average of wind stress
531      DO ij = 1, jjp1 - 2
532        zonavg = SUM(txu(ij * iip1 + 1:ij * iip1 + iim)) / iim
533        txu(ij * iip1 + 1:(ij + 1) * iip1) = zonavg
534        zonavg = SUM(tyu(ij * iip1 + 1:ij * iip1 + iim)) / iim
535        tyu(ij * iip1 + 1:(ij + 1) * iip1) = zonavg
536      END DO
537    END IF
538
539    ! Divide taux,y by f or eps, and convert to 2D u,v grids
540    ! (Arakawa C grid)
541    ! Meridional flux
542    CALL gr_scal_v(1, txu, txv) ! wind stress at v points
543    fluxm = -txv * unsfv ! in kg.s-1.m-1 (zonal distance)
544    IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
545      tcurl = (txu(1:ip1jm) - txu(iip2:ip1jmp1)) / cv(:) ! dtx/dy
546      !poles curl = 0
547      tcurl(1:iip1) = 0.
548      tcurl(ip1jmi1 + 1:ip1jm) = 0.
549      fluxm = fluxm - tcurl * unsbv
550    ENDIF
551    IF (slab_tyeq) THEN ! meridional wind forcing at equator
552      CALL gr_scal_v(1, tyu, tyv)
553      fluxm = fluxm + tyv * unsev ! in kg.s-1.m-1 (zonal distance)
554    ENDIF
555    !  apply continent mask, multiply by horiz grid dimension
556    fluxm = fluxm * cv * cuvsurcv * zmasqv
557
558    ! Zonal flux
559    IF (ekman_zonadv) THEN
560      CALL gr_scal_u(1, txu, txu) ! wind stress at u points
561      CALL gr_scal_u(1, tyu, tyu)
562      fluxz = tyu * unsfu + txu * unseu
563      !  apply continent mask, multiply by horiz grid dimension
564      fluxz = fluxz * cu * cvusurcu * zmasqu
565    END IF
566
567    !  Vertical mass flux from mass budget (divergence of horiz fluxes)
568    IF (ekman_zonadv) THEN
569      DO ij = iip2, ip1jm
570        fluxv_ek(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
571      ENDDO
572    ELSE
573      DO ij = iip2, ip1jm
574        fluxv_ek(ij) = -fluxm(ij) + fluxm(ij - iip1)
575      ENDDO
576    END IF
577    DO ij = iip1, ip1jmi1, iip1
578      fluxv_ek(ij + 1) = fluxv_ek(ij + iip1)
579    END DO
580    !  vertical mass flux at Poles
581    fluxv_ek(1) = -SUM(fluxm(1:iim))
582    fluxv_ek(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
583
584    ! Meridional heat fluxes
585    DO ij = 1, ip1jm
586      ! centered scheme
587      fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
588      fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
589    END DO
590
591    ! Zonal heat fluxes
592    ! Schema upstream
593    IF (ekman_zonadv) THEN
594      DO ij = iip2, ip1jm
595        IF (fluxz(ij)>=0.) THEN
596          fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
597          fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
598        ELSE
599          fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
600          fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
601        ENDIF
602      ENDDO
603      DO ij = iip1 * 2, ip1jmp1, iip1
604        fluxtz(ij, :) = fluxtz(ij - iim, :)
605      END DO
606    ELSE
607      fluxtz(:, :) = 0.
608    ENDIF
609
610    ! Temperature tendency, horizontal advection:
611    DO ij = iip2, ip1jm
612      dt_ek(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
613              + fluxtm(ij, :) - fluxtm(ij - iip1, :)
614    END DO
615    ! Poles
616    dt_ek(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
617    dt_ek(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)
618
619    ! ------------------------------------
620    ! GM mass fluxes and Temp tendency
621    ! ------------------------------------
622    IF (slab_gm) THEN
623      ! Vertical Temperature difference T1-T2 on v-grid points
624      CALL gr_scal_v(1, ts(:, 1) - ts(:, 2), dtz)
625      dtz(:) = MAX(dtz(:), 0.25)
626      ! Horizontal Temperature differences
627      CALL grad(1, (ts(:, 1) + ts(:, 2)) / 2., dtx, dty)
628      ! Meridional flux = -k.s (s=dyT/dzT)
629      ! Continent mask, multiply by dz/dy
630      fluxm = dty / dtz * 500. * cuvsurcv * zmasqv
631      ! slope limitation, multiply by kappa
632      fluxm = -gmkappa * SIGN(MIN(ABS(fluxm), gm_smax * cv * cuvsurcv), dty)
633      ! convert to kg/s
634      fluxm(:) = fluxm(:) * sea_den
635      ! Zonal flux = 0. (temporary)
636      fluxz(:) = 0.
637      !  Vertical mass flux from mass budget (divergence of horiz fluxes)
638      DO ij = iip2, ip1jm
639        fluxv_gm(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
640      ENDDO
641      DO ij = iip1, ip1jmi1, iip1
642        fluxv_gm(ij + 1) = fluxv_gm(ij + iip1)
643      END DO
644      !  vertical mass flux at Poles
645      fluxv_gm(1) = -SUM(fluxm(1:iim))
646      fluxv_gm(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
647
648      ! Meridional heat fluxes
649      DO ij = 1, ip1jm
650        ! centered scheme
651        fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
652        fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
653      END DO
654
655      ! Zonal heat fluxes
656      ! Schema upstream
657      DO ij = iip2, ip1jm
658        IF (fluxz(ij)>=0.) THEN
659          fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
660          fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
661        ELSE
662          fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
663          fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
664        ENDIF
665      ENDDO
666      DO ij = iip1 * 2, ip1jmp1, iip1
667        fluxtz(ij, :) = fluxtz(ij - iim, :)
668      END DO
669
670      ! Temperature tendency :
671      ! divergence of horizontal heat fluxes
672      DO ij = iip2, ip1jm
673        dt_gm(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
674                + fluxtm(ij, :) - fluxtm(ij - iip1, :)
675      END DO
676      ! Poles
677      dt_gm(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
678      dt_gm(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)
679    ELSE
680      dt_gm(:, :) = 0.
681      fluxv_gm(:) = 0.
682    ENDIF ! slab_gm
683
684    ! ------------------------------------
685    ! Temp tendency from vertical advection
686    ! Divide by cell area
687    ! ------------------------------------
688    ! vertical heat flux = mass flux * T, upstream scheme
689    DO ij = iip2, ip1jm
690      fluxv = fluxv_ek(ij) + fluxv_gm(ij) ! net flux, needed for upstream scheme
691      IF (fluxv>0.) THEN
692        dt_ek(ij, 1) = dt_ek(ij, 1) + fluxv_ek(ij) * ts(ij, 2)
693        dt_ek(ij, 2) = dt_ek(ij, 2) - fluxv_ek(ij) * ts(ij, 2)
694        dt_gm(ij, 1) = dt_gm(ij, 1) + fluxv_gm(ij) * ts(ij, 2)
695        dt_gm(ij, 2) = dt_gm(ij, 2) - fluxv_gm(ij) * ts(ij, 2)
696      ELSE
697        dt_ek(ij, 1) = dt_ek(ij, 1) + fluxv_ek(ij) * ts(ij, 1)
698        dt_ek(ij, 2) = dt_ek(ij, 2) - fluxv_ek(ij) * ts(ij, 1)
699        dt_gm(ij, 1) = dt_gm(ij, 1) + fluxv_gm(ij) * ts(ij, 1)
700        dt_gm(ij, 2) = dt_gm(ij, 2) - fluxv_gm(ij) * ts(ij, 1)
701      ENDIF
702      ! divide by cell area
703      dt_ek(ij, :) = dt_ek(ij, :) / aire(ij)
704      dt_gm(ij, :) = dt_gm(ij, :) / aire(ij)
705    END DO
706    ! North Pole
707    fluxv = fluxv_ek(1) + fluxv_gm(1)
708    IF (fluxv>0.) THEN
709      dt_ek(1, 1) = dt_ek(1, 1) + fluxv_ek(1) * ts(1, 2)
710      dt_ek(1, 2) = dt_ek(1, 2) - fluxv_ek(1) * ts(1, 2)
711      dt_gm(1, 1) = dt_gm(1, 1) + fluxv_gm(1) * ts(1, 2)
712      dt_gm(1, 2) = dt_gm(1, 2) - fluxv_gm(1) * ts(1, 2)
713    ELSE
714      dt_ek(1, 1) = dt_ek(1, 1) + fluxv_ek(1) * ts(1, 1)
715      dt_ek(1, 2) = dt_ek(1, 2) - fluxv_ek(1) * ts(1, 1)
716      dt_gm(1, 1) = dt_gm(1, 1) + fluxv_gm(1) * ts(1, 1)
717      dt_gm(1, 2) = dt_gm(1, 2) - fluxv_gm(1) * ts(1, 1)
718    ENDIF
719    dt_ek(1, :) = dt_ek(1, :) / apoln
720    dt_gm(1, :) = dt_gm(1, :) / apoln
721    ! South pole
722    fluxv = fluxv_ek(ip1jmp1) + fluxv_gm(ip1jmp1)
723    IF (fluxv>0.) THEN
724      dt_ek(ip1jmp1, 1) = dt_ek(ip1jmp1, 1) + fluxv_ek(ip1jmp1) * ts(ip1jmp1, 2)
725      dt_ek(ip1jmp1, 2) = dt_ek(ip1jmp1, 2) - fluxv_ek(ip1jmp1) * ts(ip1jmp1, 2)
726      dt_gm(ip1jmp1, 1) = dt_gm(ip1jmp1, 1) + fluxv_gm(ip1jmp1) * ts(ip1jmp1, 2)
727      dt_gm(ip1jmp1, 2) = dt_gm(ip1jmp1, 2) - fluxv_gm(ip1jmp1) * ts(ip1jmp1, 2)
728    ELSE
729      dt_ek(ip1jmp1, 1) = dt_ek(ip1jmp1, 1) + fluxv_ek(ip1jmp1) * ts(ip1jmp1, 1)
730      dt_ek(ip1jmp1, 2) = dt_ek(ip1jmp1, 2) - fluxv_ek(ip1jmp1) * ts(ip1jmp1, 1)
731      dt_gm(ip1jmp1, 1) = dt_gm(ip1jmp1, 1) + fluxv_gm(ip1jmp1) * ts(ip1jmp1, 1)
732      dt_gm(ip1jmp1, 2) = dt_gm(ip1jmp1, 2) - fluxv_gm(ip1jmp1) * ts(ip1jmp1, 1)
733    ENDIF
734    dt_ek(ip1jmp1, :) = dt_ek(ip1jmp1, :) / apols
735    dt_gm(ip1jmp1, :) = dt_gm(ip1jmp1, :) / apols
736
737    dt_ek(2:iip1, 1) = dt_ek(1, 1)
738    dt_ek(2:iip1, 2) = dt_ek(1, 2)
739    dt_gm(2:iip1, 1) = dt_gm(1, 1)
740    dt_gm(2:iip1, 2) = dt_gm(1, 2)
741    dt_ek(ip1jm + 1:ip1jmp1, 1) = dt_ek(ip1jmp1, 1)
742    dt_ek(ip1jm + 1:ip1jmp1, 2) = dt_ek(ip1jmp1, 2)
743    dt_gm(ip1jm + 1:ip1jmp1, 1) = dt_gm(ip1jmp1, 1)
744    dt_gm(ip1jm + 1:ip1jmp1, 2) = dt_gm(ip1jmp1, 2)
745
746    DO ij = iip1, ip1jmi1, iip1
747      dt_gm(ij + 1, :) = dt_gm(ij + iip1, :)
748      dt_ek(ij + 1, :) = dt_ek(ij + iip1, :)
749    END DO
750
751    ! T tendency back to 1D grid...
752    CALL gr_dyn_fi(2, iip1, jjp1, dt_ek, dt_phy_ek)
753    CALL gr_dyn_fi(2, iip1, jjp1, dt_gm, dt_phy_gm)
754
755  END SUBROUTINE slab_ekman2
756
757  SUBROUTINE slab_gmdiff(ts_phy, dt_phy)
758    ! Temperature tendency for 2-layers slab ocean
759    ! Due to Gent-McWilliams type eddy-induced advection
760
761    IMPLICIT NONE
762
763    ! Here, temperature and flux variables are on 2 layers
764    INTEGER ij
765    ! Temperature gradient, v-points
766    REAL dty((nbp_lon + 1) * (nbp_lat - 1)), dtx((nbp_lon + 1) * nbp_lat)
767    ! Vertical temperature difference, V-points
768    REAL dtz((nbp_lon + 1) * (nbp_lat - 1))
769    ! slab temperature on  1D, 2D grid
770    REAL ts_phy(klon_glo, 2), ts((nbp_lon + 1) * nbp_lat, 2)
771    ! zonal and meridional mass fluxes at u, v points (2D grid)
772    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
773    ! vertical mass flux between the 2 layers
774    REAL fluxv((nbp_lon + 1) * nbp_lat)
775    ! zonal and meridional heat fluxes
776    REAL fluxtz((nbp_lon + 1) * nbp_lat, 2)
777    REAL fluxtm((nbp_lon + 1) * (nbp_lat - 1), 2)
778    ! temperature tendency (in K.s-1.kg.m-2)
779    REAL dt((nbp_lon + 1) * nbp_lat, 2), dt_phy(klon_glo, 2)
780
781    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1
782
783    ! Grid definitions
784    iim = nbp_lon
785    iip1 = nbp_lon + 1
786    iip2 = nbp_lon + 2
787    jjp1 = nbp_lat
788    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
789    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
790    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
791
792    ! Convert temperature to 2D grid
793    CALL gr_fi_dyn(2, iip1, jjp1, ts_phy, ts)
794    ! Vertical Temperature difference T1-T2 on v-grid points
795    CALL gr_scal_v(1, ts(:, 1) - ts(:, 2), dtz)
796    dtz(:) = MAX(dtz(:), 0.25)
797    ! Horizontal Temperature differences
798    CALL grad(1, (ts(:, 1) + ts(:, 2)) / 2., dtx, dty)
799    ! Meridional flux = -k.s (s=dyT/dzT)
800    ! Continent mask, multiply by dz/dy
801    fluxm = dty / dtz * 500. * cuvsurcv * zmasqv
802    ! slope limitation, multiply by kappa
803    fluxm = -gmkappa * SIGN(MIN(ABS(fluxm), gm_smax * cv * cuvsurcv), dty)
804    ! Zonal flux = 0. (temporary)
805    fluxz(:) = 0.
806    !  Vertical mass flux from mass budget (divergence of horiz fluxes)
807    DO ij = iip2, ip1jm
808      fluxv(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
809    ENDDO
810    DO ij = iip1, ip1jmi1, iip1
811      fluxv(ij + 1) = fluxv(ij + iip1)
812    END DO
813    !  vertical mass flux at Poles
814    fluxv(1) = -SUM(fluxm(1:iim))
815    fluxv(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
816    fluxv = fluxv
817
818    ! Meridional heat fluxes
819    DO ij = 1, ip1jm
820      ! centered scheme
821      fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
822      fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
823    END DO
824
825    ! Zonal heat fluxes
826    ! Schema upstream
827    DO ij = iip2, ip1jm
828      IF (fluxz(ij)>=0.) THEN
829        fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
830        fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
831      ELSE
832        fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
833        fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
834      ENDIF
835    ENDDO
836    DO ij = iip1 * 2, ip1jmp1, iip1
837      fluxtz(ij, :) = fluxtz(ij - iim, :)
838    END DO
839
840    ! Temperature tendency :
841    DO ij = iip2, ip1jm
842      ! divergence of horizontal heat fluxes
843      dt(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
844              + fluxtm(ij, :) - fluxtm(ij - iip1, :)
845      ! + vertical heat flux (mass flux * T, upstream scheme)
846      IF (fluxv(ij)>0.) THEN
847        dt(ij, 1) = dt(ij, 1) + fluxv(ij) * ts(ij, 2)
848        dt(ij, 2) = dt(ij, 2) - fluxv(ij) * ts(ij, 2)
849      ELSE
850        dt(ij, 1) = dt(ij, 1) + fluxv(ij) * ts(ij, 1)
851        dt(ij, 2) = dt(ij, 2) - fluxv(ij) * ts(ij, 1)
852      ENDIF
853      ! divide by cell area
854      dt(ij, :) = dt(ij, :) / aire(ij)
855    END DO
856    DO ij = iip1, ip1jmi1, iip1
857      dt(ij + 1, :) = dt(ij + iip1, :)
858    END DO
859    ! Poles
860    dt(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
861    IF (fluxv(1)>0.) THEN
862      dt(1, 1) = dt(1, 1) + fluxv(1) * ts(1, 2)
863      dt(1, 2) = dt(1, 2) - fluxv(1) * ts(1, 2)
864    ELSE
865      dt(1, 1) = dt(1, 1) + fluxv(1) * ts(1, 1)
866      dt(1, 2) = dt(1, 2) - fluxv(1) * ts(1, 1)
867    ENDIF
868    dt(1, :) = dt(1, :) / apoln
869    dt(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)
870    IF (fluxv(ip1jmp1)>0.) THEN
871      dt(ip1jmp1, 1) = dt(ip1jmp1, 1) + fluxv(ip1jmp1) * ts(ip1jmp1, 2)
872      dt(ip1jmp1, 2) = dt(ip1jmp1, 2) - fluxv(ip1jmp1) * ts(ip1jmp1, 2)
873    ELSE
874      dt(ip1jmp1, 1) = dt(ip1jmp1, 1) + fluxv(ip1jmp1) * ts(ip1jmp1, 1)
875      dt(ip1jmp1, 2) = dt(ip1jmp1, 2) - fluxv(ip1jmp1) * ts(ip1jmp1, 1)
876    ENDIF
877    dt(ip1jmp1, :) = dt(ip1jmp1, :) / apols
878    dt(2:iip1, 1) = dt(1, 1)
879    dt(2:iip1, 2) = dt(1, 2)
880    dt(ip1jm + 1:ip1jmp1, 1) = dt(ip1jmp1, 1)
881    dt(ip1jm + 1:ip1jmp1, 2) = dt(ip1jmp1, 2)
882
883    ! T tendency back to 1D grid...
884    CALL gr_dyn_fi(2, iip1, jjp1, dt, dt_phy)
885
886  END SUBROUTINE slab_gmdiff
887
888  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
889
890  SUBROUTINE gr_fi_dyn(nfield, im, jm, pfi, pdyn)
891    ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
892    USE lmdz_ssum_scopy, ONLY: scopy
893
894    IMPLICIT NONE
895
896    INTEGER, INTENT(IN) :: im, jm, nfield
897    REAL, INTENT(IN) :: pfi(klon_glo, nfield) ! on 1D grid
898    REAL, INTENT(OUT) :: pdyn(im, jm, nfield) ! on 2D grid
899
900    INTEGER :: i, j, ifield, ig
901
902    DO ifield = 1, nfield
903      ! Handle poles
904      DO i = 1, im
905        pdyn(i, 1, ifield) = pfi(1, ifield)
906        pdyn(i, jm, ifield) = pfi(klon_glo, ifield)
907      ENDDO
908      ! Other points
909      DO j = 2, jm - 1
910        ig = 2 + (j - 2) * (im - 1)
911        CALL SCOPY(im - 1, pfi(ig, ifield), 1, pdyn(1, j, ifield), 1)
912        pdyn(im, j, ifield) = pdyn(1, j, ifield)
913      ENDDO
914    ENDDO ! of DO ifield=1,nfield
915
916  END SUBROUTINE gr_fi_dyn
917
918  SUBROUTINE gr_dyn_fi(nfield, im, jm, pdyn, pfi)
919    ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
920    USE lmdz_ssum_scopy, ONLY: scopy
921    IMPLICIT NONE
922
923    INTEGER, INTENT(IN) :: im, jm, nfield
924    REAL, INTENT(IN) :: pdyn(im, jm, nfield) ! on 2D grid
925    REAL, INTENT(OUT) :: pfi(klon_glo, nfield) ! on 1D grid
926
927    INTEGER j, ifield, ig
928
929    CHARACTER (len = 20) :: modname = 'slab_heat_transp'
930    CHARACTER (len = 80) :: abort_message
931
932    ! Sanity check:
933    IF(klon_glo/=2 + (jm - 2) * (im - 1)) THEN
934      abort_message = "gr_dyn_fi error, wrong sizes"
935      CALL abort_physic(modname, abort_message, 1)
936    ENDIF
937
938    ! Handle poles
939    CALL SCOPY(nfield, pdyn, im * jm, pfi, klon_glo)
940    CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(klon_glo, 1), klon_glo)
941    ! Other points
942    DO ifield = 1, nfield
943      DO j = 2, jm - 1
944        ig = 2 + (j - 2) * (im - 1)
945        CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1)
946      ENDDO
947    ENDDO
948
949  END SUBROUTINE gr_dyn_fi
950
951  SUBROUTINE  grad(klevel, pg, pgx, pgy)
952    ! compute the covariant components pgx,pgy of the gradient of pg
953    ! pgx = d(pg)/dx * delta(x) = delta(pg)
954    IMPLICIT NONE
955
956    INTEGER, INTENT(IN) :: klevel
957    REAL, INTENT(IN) :: pg((nbp_lon + 1) * nbp_lat, klevel)
958    REAL, INTENT(OUT) :: pgx((nbp_lon + 1) * nbp_lat, klevel)
959    REAL, INTENT(OUT) :: pgy((nbp_lon + 1) * (nbp_lat - 1), klevel)
960
961    INTEGER :: l, ij
962    INTEGER :: iim, iip1, ip1jm, ip1jmp1
963
964    iim = nbp_lon
965    iip1 = nbp_lon + 1
966    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
967    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
968
969    DO l = 1, klevel
970      DO ij = 1, ip1jmp1 - 1
971        pgx(ij, l) = pg(ij + 1, l) - pg(ij, l)
972      ENDDO
973      ! correction for pgx(ip1,j,l) ...
974      ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
975      DO ij = iip1, ip1jmp1, iip1
976        pgx(ij, l) = pgx(ij - iim, l)
977      ENDDO
978      DO ij = 1, ip1jm
979        pgy(ij, l) = pg(ij, l) - pg(ij + iip1, l)
980      ENDDO
981    ENDDO
982
983  END SUBROUTINE grad
984
985  SUBROUTINE diverg(klevel, x, y, div)
986    ! computes the divergence of a vector field of components
987    ! x,y. x and y being covariant components
988    USE lmdz_ssum_scopy, ONLY: ssum
989
990    IMPLICIT NONE
991
992    INTEGER, INTENT(IN) :: klevel
993    REAL, INTENT(IN) :: x((nbp_lon + 1) * nbp_lat, klevel)
994    REAL, INTENT(IN) :: y((nbp_lon + 1) * (nbp_lat - 1), klevel)
995    REAL, INTENT(OUT) :: div((nbp_lon + 1) * nbp_lat, klevel)
996
997    INTEGER :: l, ij
998    INTEGER :: iim, iip1, iip2, ip1jm, ip1jmp1, ip1jmi1
999
1000    REAL :: aiy1(nbp_lon + 1), aiy2(nbp_lon + 1)
1001    REAL :: sumypn, sumyps
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    USE lmdz_ssum_scopy, ONLY: scopy
1107
1108    IMPLICIT NONE
1109
1110    INTEGER, INTENT(IN) :: nx
1111    REAL, INTENT(OUT) :: x_u((nbp_lon + 1) * nbp_lat, nx)
1112    REAL, INTENT(IN) :: x_scal((nbp_lon + 1) * nbp_lat, nx)
1113
1114    INTEGER :: l, ij
1115    INTEGER :: iip1, jjp1, ip1jmp1
1116
1117    iip1 = nbp_lon + 1
1118    jjp1 = nbp_lat
1119    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
1120
1121    DO l = 1, nx
1122      DO ij = 1, ip1jmp1 - 1
1123        x_u(ij, l) = &
1124                (aire(ij) * x_scal(ij, l) + aire(ij + 1) * x_scal(ij + 1, l)) &
1125                        / (aire(ij) + aire(ij + 1))
1126      ENDDO
1127    ENDDO
1128
1129    CALL SCOPY(nx * jjp1, x_u(1, 1), iip1, x_u(iip1, 1), iip1)
1130
1131  END SUBROUTINE gr_scal_u
1132
1133END MODULE slab_heat_transp_mod
Note: See TracBrowser for help on using the repository browser.