source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/slab_heat_transp_mod.F90 @ 3603

Last change on this file since 3603 was 3406, checked in by jghattas, 6 years ago

Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
File size: 36.2 KB
Line 
1!
2MODULE slab_heat_transp_mod
3!
4! Slab ocean : temperature tendencies due to horizontal diffusion
5! and / or Ekman transport
6
7USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo
8IMPLICIT 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)
86    !!USE comconst_mod, ONLY: omeg, rad
87    ! number of points in lon, lat
88    IMPLICIT NONE
89    ! Routine copies some geometry variables from the dynamical core
90    ! see global vars for meaning
91    INTEGER,INTENT(IN) :: ip1jm
92    INTEGER,INTENT(IN) :: ip1jmp1
93    REAL,INTENT(IN) :: unsairez_(ip1jm)
94    REAL,INTENT(IN) :: fext_(ip1jm)
95    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
96    REAL,INTENT(IN) :: cu_(ip1jmp1)
97    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
98    REAL,INTENT(IN) :: cv_(ip1jm)
99    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
100    REAL,INTENT(IN) :: aire_(ip1jmp1)
101    REAL,INTENT(IN) :: apoln_
102    REAL,INTENT(IN) :: apols_
103    REAL,INTENT(IN) :: aireu_(ip1jmp1)
104    REAL,INTENT(IN) :: airev_(ip1jm)
105    REAL,INTENT(IN) :: rlatv(nbp_lat-1)
106    REAL            :: omeg = 7.272205e-05
107    REAL            :: rad = 6371229.
108
109    ! Sanity check on dimensions
110    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
111        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
112      write(*,*) "ini_slab_transp_geom Error: wrong array sizes"
113      stop
114    endif
115! Allocations could be done only on master process/thread...
116    allocate(unsairez(ip1jm))
117    unsairez(:)=unsairez_(:)
118    allocate(fext(ip1jm))
119    fext(:)=fext_(:)
120    allocate(unsaire(ip1jmp1))
121    unsaire(:)=unsaire_(:)
122    allocate(cu(ip1jmp1))
123    cu(:)=cu_(:)
124    allocate(cuvsurcv(ip1jm))
125    cuvsurcv(:)=cuvsurcv_(:)
126    allocate(cv(ip1jm))
127    cv(:)=cv_(:)
128    allocate(cvusurcu(ip1jmp1))
129    cvusurcu(:)=cvusurcu_(:)
130    allocate(aire(ip1jmp1))
131    aire(:)=aire_(:)
132    apoln=apoln_
133    apols=apols_
134    allocate(aireu(ip1jmp1))
135    aireu(:)=aireu_(:)
136    allocate(airev(ip1jm))
137    airev(:)=airev_(:)
138    allocate(beta(nbp_lat-1))
139    beta(:)=2*omeg*cos(rlatv(:))/rad
140
141  END SUBROUTINE ini_slab_transp_geom
142
143  SUBROUTINE ini_slab_transp(zmasq)
144
145!    USE ioipsl_getin_p_mod, only: getin_p
146    USE IOIPSL, ONLY : getin
147    IMPLICIT NONE
148
149    REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent
150    REAL zmasq_2d((nbp_lon+1)*nbp_lat)
151    REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter
152    REAL eps ! epsilon friction timescale (s-1)
153    INTEGER :: slab_ekman
154    INTEGER i
155    INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
156
157! Some definition for grid size
158    ip1jm=(nbp_lon+1)*(nbp_lat-1)
159    ip1jmp1=(nbp_lon+1)*nbp_lat
160    iim=nbp_lon
161    iip1=nbp_lon+1
162    jjp1=nbp_lat
163    ip1jm=(nbp_lon+1)*(nbp_lat-1)
164    ip1jmp1=(nbp_lon+1)*nbp_lat
165
166! Options for Heat transport
167    ! Alpha variable?
168      alpha_var=.FALSE.
169      CALL getin('slab_alphav',alpha_var)
170      print *,'alpha variable',alpha_var
171!  centered ou upstream scheme for meridional transport
172      slab_upstream=.FALSE.
173      CALL getin('slab_upstream',slab_upstream)
174      print *,'upstream slab scheme',slab_upstream
175! Sverdrup balance at equator ?
176      slab_sverdrup=.TRUE.
177      CALL getin('slab_sverdrup',slab_sverdrup)
178      print *,'Sverdrup balance',slab_sverdrup
179! Use tauy for meridional flux at equator ?
180      slab_tyeq=.TRUE.
181      CALL getin('slab_tyeq',slab_tyeq)
182      print *,'Tauy forcing at equator',slab_tyeq
183! Use tauy for meridional flux at equator ?
184      ekman_zonadv=.TRUE.
185      CALL getin('slab_ekman_zonadv',ekman_zonadv)
186      print *,'Use Ekman flow in zonal direction',ekman_zonadv
187! Use tauy for meridional flux at equator ?
188      ekman_zonavg=.FALSE.
189      CALL getin('slab_ekman_zonavg',ekman_zonavg)
190      print *,'Use zonally-averaged wind stress ?',ekman_zonavg
191! Value of alpha
192      alpham=2./3.
193      CALL getin('slab_alpha',alpham)
194      print *,'slab_alpha',alpham
195! GM k coefficient (m2/s) for 2-layers
196      gmkappa=1000.
197      CALL getin('slab_gmkappa',gmkappa)
198      print *,'slab_gmkappa',gmkappa
199! GM k coefficient (m2/s) for 2-layers
200      gm_smax=2e-3
201      CALL getin('slab_gm_smax',gm_smax)
202      print *,'slab_gm_smax',gm_smax
203! -----------------------------------------------------------
204! Define ocean / continent mask (no flux into continent cell)
205    allocate(zmasqu(ip1jmp1))
206    allocate(zmasqv(ip1jm))
207    zmasqu=1.
208    zmasqv=1.
209
210    ! convert mask to 2D grid
211    CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d)
212    ! put flux mask to 0 at boundaries of continent cells
213    DO i=1,ip1jmp1-1
214      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
215              zmasqu(i)=0.
216      ENDIF
217    END DO
218    DO i=iip1,ip1jmp1,iip1
219      zmasqu(i)=zmasqu(i-iim)
220    END DO
221    DO i=1,ip1jm
222      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN
223              zmasqv(i)=0.
224      END IF
225    END DO
226
227! -----------------------------------------------------------
228! Coriolis and friction for Ekman transport
229    slab_ekman=2
230    CALL getin("slab_ekman",slab_ekman)
231    IF (slab_ekman.GT.0) THEN
232      allocate(unsfv(ip1jm))
233      allocate(unsev(ip1jm))
234      allocate(unsfu(ip1jmp1))
235      allocate(unseu(ip1jmp1))
236      allocate(unsbv(ip1jm))
237
238      eps=1e-5 ! Drag
239      CALL getin('slab_eps',eps)
240      print *,'epsilon=',eps
241      ff=fext*unsairez ! Coriolis
242      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
243      ! on 2D grid v points
244      ! Compute correction factor [0 1] near the equator (f<<eps)
245      IF (slab_sverdrup) THEN
246         ! New formulation, sharper near equator, when eps gives Rossby Radius
247         DO i=1,ip1jm
248           unsev(i)=exp(-ff(i)*ff(i)/eps**2)
249         ENDDO
250      ELSE
251         DO i=1,ip1jm
252           unsev(i)=eps**2/(ff(i)*ff(i)+eps**2)
253         ENDDO
254      END IF ! slab_sverdrup
255      ! 1/beta
256      DO i=1,jjp1-1
257        unsbv((i-1)*iip1+1:i*iip1)=unsev((i-1)*iip1+1:i*iip1)/beta(i)
258      END DO
259      ! 1/f
260      ff=SIGN(MAX(ABS(ff),eps/100.),ff) ! avoid value 0 at equator...
261      DO i=1,ip1jm
262        unsfv(i)=(1.-unsev(i))/ff(i)
263      END DO
264      ! compute values on 2D u grid
265      ! 1/eps
266      unsev(:)=unsev(:)/eps
267      CALL gr_v_scal(1,unsfv,unsfu)
268      CALL gr_v_scal(1,unsev,unseu)
269    END IF
270 
271  END SUBROUTINE ini_slab_transp
272
273  SUBROUTINE divgrad_phy(nlevs,temp,delta)
274! Computes temperature tendency due to horizontal diffusion :
275! T Laplacian, later multiplied by diffusion coef and time-step
276
277    IMPLICIT NONE
278 
279    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
280    REAL, INTENT(IN)   :: temp(klon_glo,nlevs) ! slab temperature
281    REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.)
282    REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs)
283    REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
284    INTEGER :: ll,iip1,jjp1
285 
286    iip1=nbp_lon+1
287    jjp1=nbp_lat
288   
289    ! transpose temp to 2D horiz. grid
290    CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d)
291    ! computes gradient (proportional to heat flx)
292    CALL grad(nlevs,delta_2d,ghx,ghy)
293    ! put flux to 0 at ocean / continent boundary
294    DO ll=1,nlevs
295        ghx(:,ll)=ghx(:,ll)*zmasqu
296        ghy(:,ll)=ghy(:,ll)*zmasqv
297    END DO
298    ! flux divergence
299    CALL diverg(nlevs,ghx,ghy,delta_2d)
300    ! laplacian back to 1D grid
301    CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta)
302 
303    RETURN
304  END SUBROUTINE divgrad_phy
305
306  SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy)
307! 1.5 Layer Ekman transport temperature tendency
308! note : tendency dt later multiplied by (delta t)/(rho.H)
309! to convert from divergence of heat fluxes to T
310
311      IMPLICIT NONE
312     
313      ! tx, ty : wind stress (different grids)
314      ! fluxm, fluz : mass *or heat* fluxes
315      ! dt : temperature tendency
316      INTEGER ij
317
318      ! ts surface temp, td deep temp (diagnosed)
319      REAL ts_phy(klon_glo)
320      REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat)
321      ! zonal and meridional wind stress
322      REAL tx_phy(klon_glo),ty_phy(klon_glo)
323      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
324      REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1))
325      REAL tcurl((nbp_lon+1)*(nbp_lat-1))
326      ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
327      REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1))
328      ! vertical  and absolute mass fluxes (to estimate alpha)
329      REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1))
330      ! temperature tendency
331      REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo)
332      REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef
333
334      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
335
336! Grid definitions
337      iim=nbp_lon
338      iip1=nbp_lon+1
339      iip2=nbp_lon+2
340      jjp1=nbp_lat
341      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
342      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
343      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
344
345! Convert taux,y to 2D  scalar grid
346! Note: 2D grid size = iim*jjm. iip1=iim+1
347! First and last points in zonal direction are the same
348! we use 1 index ij from 1 to (iim+1)*(jjm+1)
349      ! north and south poles
350      tx_phy(1)=0.
351      tx_phy(klon_glo)=0.
352      ty_phy(1)=0.
353      ty_phy(klon_glo)=0.
354      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
355      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
356! convert to u,v grid (Arakawa C)
357! Multiply by f or eps to get mass flux
358      ! Meridional mass flux
359      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
360      IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
361        tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:)
362        fluxm=-tcurl*unsbv-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
363      ELSE
364        CALL gr_scal_v(1,tyu,tyv)
365        fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
366      ENDIF
367      ! Zonal mass flux
368      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
369      CALL gr_scal_u(1,tyu,tyu)
370      fluxz=tyu*unsfu+txu*unseu
371           
372! Correct flux: continent mask and horiz grid size
373      ! multiply m-flux by mask and dx: flux in kg.s-1
374      fluxm=fluxm*cv*cuvsurcv*zmasqv
375      ! multiply z-flux by mask and dy: flux in kg.s-1
376      fluxz=fluxz*cu*cvusurcu*zmasqu
377
378! Compute vertical  and absolute mass flux (for variable alpha)
379      IF (alpha_var) THEN
380        DO ij=iip2,ip1jm
381        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
382        fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) &
383               +ABS(fluxm(ij))+ABS(fluxm(ij-iip1))
384        ENDDO
385        DO ij=iip1,ip1jmi1,iip1
386            fluxt(ij+1)=fluxt(ij+iip1)
387            fluxv(ij+1)=fluxv(ij+iip1)
388        END DO
389        fluxt(1)=SUM(ABS(fluxm(1:iim)))
390        fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1)))
391        fluxv(1)=-SUM(fluxm(1:iim))
392        fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
393        fluxt=MAX(fluxt,1.e-10)
394      ENDIF
395
396! Compute alpha coefficient.
397! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
398      IF (alpha_var) THEN
399          ! increase alpha (and Tdeep) in downwelling regions
400          ! and decrease in upwelling regions
401          ! to avoid "hot spots" where there is surface convergence
402          DO ij=iip2,ip1jm
403              alpha(ij)=alpham-fluxv(ij)/fluxt(ij)*(1.-alpham)
404          ENDDO
405          alpha(1:iip1)=alpham-fluxv(1)/fluxt(1)*(1.-alpham)
406          alpha(ip1jm+1:ip1jmp1)=alpham-fluxv(ip1jmp1)/fluxt(ip1jmp1)*(1.-alpham)
407      ELSE
408          alpha(:)=alpham
409          ! Tsurf-Tdeep ~ 10° in the Tropics
410      ENDIF
411
412! Estimate deep temperature
413      CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts)
414      DO ij=1,ip1jmp1
415         td(ij)=271.35+(ts(ij)-271.35)*alpha(ij)
416         td(ij)=MIN(td(ij),ts(ij))
417      END DO
418       
419! Meridional heat flux: multiply mass flux by (ts-td)
420! flux in K.kg.s-1
421      IF (slab_upstream) THEN
422        ! upstream scheme to avoid hot spots
423        DO ij=1,ip1jm
424        IF (fluxm(ij).GE.0.) THEN
425           fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij))
426        ELSE
427           fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1))
428        END IF
429        END DO
430      ELSE
431        ! centered scheme better in mid-latitudes
432        DO ij=1,ip1jm
433          fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2.
434        END DO
435      ENDIF
436
437! Zonal heat flux
438! upstream scheme
439      DO ij=iip2,ip1jm
440          fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2.
441      END DO
442      DO ij=iip1*2,ip1jmp1,iip1
443           fluxz(ij)=fluxz(ij-iim)
444      END DO
445
446! temperature tendency = divergence of heat fluxes
447! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
448      DO ij=iip2,ip1jm
449         dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) &
450               /aire(ij) ! aire : grid area
451      END DO
452      DO ij=iip1,ip1jmi1,iip1
453         dt(ij+1)=dt(ij+iip1)
454      END DO
455! special treatment at the Poles
456      dt(1)=SUM(fluxm(1:iim))/apoln
457      dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols
458      dt(2:iip1)=dt(1)
459      dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1)
460
461! tendencies back to 1D grid
462      CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy)
463
464      RETURN
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).GE.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).GE.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.GT.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.GT.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.GT.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      RETURN
756  END SUBROUTINE slab_ekman2
757
758  SUBROUTINE slab_gmdiff(ts_phy,dt_phy)
759! Temperature tendency for 2-layers slab ocean
760! Due to Gent-McWilliams type eddy-induced advection
761     
762      IMPLICIT NONE
763     
764      ! Here, temperature and flux variables are on 2 layers
765      INTEGER ij
766      ! Temperature gradient, v-points
767      REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat)
768      ! Vertical temperature difference, V-points
769      REAL dtz((nbp_lon+1)*(nbp_lat-1))
770      ! slab temperature on  1D, 2D grid
771      REAL ts_phy(klon_glo,2),ts((nbp_lon+1)*nbp_lat,2)
772      ! zonal and meridional mass fluxes at u, v points (2D grid)
773      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
774      ! vertical mass flux between the 2 layers
775      REAL fluxv((nbp_lon+1)*nbp_lat)
776      ! zonal and meridional heat fluxes
777      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
778      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
779      ! temperature tendency (in K.s-1.kg.m-2)
780      REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2)
781
782      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
783
784! Grid definitions
785      iim=nbp_lon
786      iip1=nbp_lon+1
787      iip2=nbp_lon+2
788      jjp1=nbp_lat
789      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
790      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
791      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
792
793! Convert temperature to 2D grid
794      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
795! Vertical Temperature difference T1-T2 on v-grid points
796      CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz)
797      dtz(:)=MAX(dtz(:),0.25)
798! Horizontal Temperature differences
799      CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty)
800! Meridional flux = -k.s (s=dyT/dzT)
801! Continent mask, multiply by dz/dy
802      fluxm=dty/dtz*500.*cuvsurcv*zmasqv
803! slope limitation, multiply by kappa
804      fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty)
805! Zonal flux = 0. (temporary)
806      fluxz(:)=0.
807!  Vertical mass flux from mass budget (divergence of horiz fluxes)
808      DO ij=iip2,ip1jm
809        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
810      ENDDO
811      DO ij=iip1,ip1jmi1,iip1
812         fluxv(ij+1)=fluxv(ij+iip1)
813      END DO
814!  vertical mass flux at Poles
815      fluxv(1)=-SUM(fluxm(1:iim))     
816      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
817      fluxv=fluxv
818
819! Meridional heat fluxes
820      DO ij=1,ip1jm
821          ! centered scheme
822          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
823          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
824      END DO
825
826! Zonal heat fluxes
827! Schema upstream     
828      DO ij=iip2,ip1jm
829      IF (fluxz(ij).GE.0.) THEN
830             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
831             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
832      ELSE
833             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
834             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
835      ENDIF
836      ENDDO
837      DO ij=iip1*2,ip1jmp1,iip1
838             fluxtz(ij,:)=fluxtz(ij-iim,:)
839      END DO
840                   
841! Temperature tendency :
842      DO ij=iip2,ip1jm
843! divergence of horizontal heat fluxes
844         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
845                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
846! + vertical heat flux (mass flux * T, upstream scheme)
847         IF (fluxv(ij).GT.0.) THEN
848           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
849           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
850         ELSE
851           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
852           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
853         ENDIF
854         ! divide by cell area
855         dt(ij,:)=dt(ij,:)/aire(ij)
856      END DO
857      DO ij=iip1,ip1jmi1,iip1
858         dt(ij+1,:)=dt(ij+iip1,:)
859      END DO
860! Poles
861      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
862        IF (fluxv(1).GT.0.) THEN
863          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
864          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
865        ELSE
866          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
867          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
868        ENDIF
869      dt(1,:)=dt(1,:)/apoln
870      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
871       IF (fluxv(ip1jmp1).GT.0.) THEN
872         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
873         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
874       ELSE
875         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
876         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
877       ENDIF
878      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
879      dt(2:iip1,1)=dt(1,1)
880      dt(2:iip1,2)=dt(1,2)
881      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
882      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
883
884! T tendency back to 1D grid...
885      CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy)
886
887      RETURN
888  END SUBROUTINE slab_gmdiff
889
890!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
891
892  SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn)
893  ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
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  IMPLICIT NONE
921
922  INTEGER,INTENT(IN) :: im,jm,nfield
923  REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid
924  REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid
925
926  INTEGER j,ifield,ig
927
928  ! Sanity check:
929  IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN
930    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
931    STOP
932  ENDIF
933
934  ! Handle poles
935  CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo)
936  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo)
937  ! Other points
938  DO ifield=1,nfield
939    DO j=2,jm-1
940      ig=2+(j-2)*(im-1)
941      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
942    ENDDO
943  ENDDO
944
945  END SUBROUTINE gr_dyn_fi
946
947  SUBROUTINE  grad(klevel,pg,pgx,pgy)
948  ! compute the covariant components pgx,pgy of the gradient of pg
949  ! pgx = d(pg)/dx * delta(x) = delta(pg)
950  IMPLICIT NONE
951
952  INTEGER,INTENT(IN) :: klevel
953  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
954  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
955  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
956
957  INTEGER :: l,ij
958  INTEGER :: iim,iip1,ip1jm,ip1jmp1
959
960  iim=nbp_lon
961  iip1=nbp_lon+1
962  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
963  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
964
965  DO l=1,klevel
966    DO ij=1,ip1jmp1-1
967      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
968    ENDDO
969    ! correction for pgx(ip1,j,l) ...
970    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
971    DO ij=iip1,ip1jmp1,iip1
972      pgx(ij,l)=pgx(ij-iim,l)
973    ENDDO
974    DO ij=1,ip1jm
975      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
976    ENDDO
977  ENDDO
978
979  END SUBROUTINE grad
980
981  SUBROUTINE diverg(klevel,x,y,div)
982  ! computes the divergence of a vector field of components
983  ! x,y. x and y being covariant components
984  IMPLICIT NONE
985
986  INTEGER,INTENT(IN) :: klevel
987  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
988  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
989  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
990
991  INTEGER :: l,ij
992  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
993
994  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
995  REAL :: sumypn,sumyps
996  REAL,EXTERNAL :: SSUM
997
998  iim=nbp_lon
999  iip1=nbp_lon+1
1000  iip2=nbp_lon+2
1001  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
1002  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
1003  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
1004
1005  DO l=1,klevel
1006    DO ij=iip2,ip1jm-1
1007      div(ij+1,l)= &
1008        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
1009        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
1010    ENDDO
1011    ! correction for div(1,j,l) ...
1012    ! ... div(1,j,l)= div(iip1,j,l) ...
1013    DO ij=iip2,ip1jm,iip1
1014      div(ij,l)=div(ij+iim,l)
1015    ENDDO
1016    ! at the poles
1017    DO ij=1,iim
1018      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
1019      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
1020    ENDDO
1021    sumypn=SSUM(iim,aiy1,1)/apoln
1022    sumyps=SSUM(iim,aiy2,1)/apols
1023    DO ij=1,iip1
1024      div(ij,l)=-sumypn
1025      div(ij+ip1jm,l)=sumyps
1026    ENDDO
1027    ! End (poles)
1028  ENDDO ! of DO l=1,klevel
1029
1030  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
1031  DO l=1,klevel
1032    DO ij=iip2,ip1jm
1033      div(ij,l)=div(ij,l)*unsaire(ij)
1034    ENDDO
1035  ENDDO
1036
1037  END SUBROUTINE diverg
1038
1039  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
1040  ! convert values from v points to scalar points on C-grid
1041  ! used to  compute unsfu, unseu (u points, but depends only on latitude)
1042  IMPLICIT NONE
1043
1044  INTEGER,INTENT(IN) :: nx ! number of levels or fields
1045  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
1046  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
1047
1048  INTEGER :: l,ij
1049  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
1050
1051  iip1=nbp_lon+1
1052  iip2=nbp_lon+2
1053  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
1054  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
1055
1056  DO l=1,nx
1057    DO ij=iip2,ip1jm
1058      x_scal(ij,l)= &
1059                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
1060                  /(airev(ij-iip1)+airev(ij))
1061    ENDDO
1062    DO ij=1,iip1
1063      x_scal(ij,l)=0.
1064    ENDDO
1065    DO ij=ip1jm+1,ip1jmp1
1066      x_scal(ij,l)=0.
1067    ENDDO
1068  ENDDO
1069
1070  END SUBROUTINE gr_v_scal
1071 
1072  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
1073  ! convert values from scalar points to v points on C-grid
1074  ! used to compute wind stress at V points
1075  IMPLICIT NONE
1076
1077  INTEGER,INTENT(IN) :: nx ! number of levels or fields
1078  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
1079  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
1080
1081  INTEGER :: l,ij
1082  INTEGER :: iip1,ip1jm
1083
1084  iip1=nbp_lon+1
1085  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
1086
1087      DO l=1,nx
1088        DO ij=1,ip1jm
1089          x_v(ij,l)= &
1090            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
1091            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
1092            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
1093        ENDDO
1094      ENDDO
1095
1096  END SUBROUTINE gr_scal_v
1097
1098  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
1099  ! convert values from scalar points to U points on C-grid
1100  ! used to compute wind stress at U points
1101  IMPLICIT NONE
1102
1103  INTEGER,INTENT(IN) :: nx
1104  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
1105  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
1106
1107  INTEGER :: l,ij
1108  INTEGER :: iip1,jjp1,ip1jmp1
1109
1110  iip1=nbp_lon+1
1111  jjp1=nbp_lat
1112  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
1113
1114  DO l=1,nx
1115     DO ij=1,ip1jmp1-1
1116        x_u(ij,l)= &
1117         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
1118         /(aire(ij)+aire(ij+1))
1119     ENDDO
1120  ENDDO
1121
1122  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
1123
1124  END SUBROUTINE gr_scal_u
1125
1126END MODULE slab_heat_transp_mod
Note: See TracBrowser for help on using the repository browser.