source: LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90 @ 5423

Last change on this file since 5423 was 3002, checked in by Ehouarn Millour, 7 years ago

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