source: trunk/LMDZ.GENERIC/libf/phystd/slab_heat_transp_mod.F90 @ 3562

Last change on this file since 3562 was 3100, checked in by bhatnags, 14 months ago

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM
File size: 31.5 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 averaged wind stress
50  !$OMP THREADPRIVATE(ekman_zonavg)
51
52  REAL,SAVE :: alpham
53  !$OMP THREADPRIVATE(alpham)
54  REAL,SAVE :: gmkappa ! GM scheme coefficient
55  !$OMP THREADPRIVATE(gmkappa)
56  REAL,SAVE :: gm_smax ! GM scheme limiter (max slope)
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 comcstfi_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    beta(:)=rlatv(:) ! Temporarily storing rlatv in beta, will be updated in ini_slab_transp
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    USE comcstfi_mod, ONLY: omeg, rad
146    IMPLICIT NONE
147
148    REAL, INTENT (IN) :: zmasq(klon_glo) ! ocean / continent mask, 1=continent
149    REAL zmasq_2d((nbp_lon+1)*nbp_lat)
150    REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter
151    REAL eps ! epsilon friction timescale (s-1)
152!    INTEGER :: slab_ekman
153    LOGICAL :: 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? (for 1.5-layer Ekman)
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 (wind stress curl) 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! Zonal Ekman transport ? (Or meridional only)
184      ekman_zonadv=.TRUE.
185      CALL getin('slab_ekman_zonadv',ekman_zonadv)
186      print *,'Use Ekman flow in zonal direction',ekman_zonadv
187! Zonally average the Ekman flux ?
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 (Ekman 1.5-layer Tdeep computation)
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 max slope 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    slab_ekman=.TRUE.
231    CALL getin("slab_ekman",slab_ekman)
232!    IF (slab_ekman.GT.0) THEN
233    IF (slab_ekman) THEN
234      allocate(unsfv(ip1jm))
235      allocate(unsev(ip1jm))
236      allocate(unsfu(ip1jmp1))
237      allocate(unseu(ip1jmp1))
238      allocate(unsbv(ip1jm))
239
240      ! Drage parameter (in s-1)
241      eps=1e-5
242      CALL getin('slab_eps',eps)
243      print *,'epsilon=',eps
244      ! Coriolis paraleter
245      ff=fext*unsairez
246
247      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
248      ! on 2D grid v points
249      IF (slab_sverdrup) THEN
250      ! unsev factor [0 1] for transition Ekman / Sverdrup near the equator (f<<eps)
251      ! f/eps gives Rossby Radius
252         DO i=1,ip1jm
253           unsev(i)=exp(-ff(i)*ff(i)/eps**2)
254         ENDDO
255      ELSE
256         ! Drag vs Coriolis balance
257         DO i=1,ip1jm
258           unsev(i)=eps**2/(ff(i)*ff(i)+eps**2)
259         ENDDO
260      END IF ! slab_sverdrup
261
262!      beta(:)=2*omeg*cos(rlatv(:))/rad
263!!!!! Here, beta contains rlatv from ini_slab_transp_geom !!!!!
264      beta(:)=2*omeg*cos(beta(:))/rad
265!!!!!
266      ! 1/beta (derivative of Coriolos parameter)
267      DO i=1,jjp1-1
268        unsbv((i-1)*iip1+1:i*iip1)=unsev((i-1)*iip1+1:i*iip1)/beta(i)
269      END DO
270      ! 1/f
271      ff=SIGN(MAX(ABS(ff),eps/100.),ff) ! avoid value 0 at equator...
272      DO i=1,ip1jm
273        unsfv(i)=(1.-unsev(i))/ff(i)
274      END DO
275      ! compute values on 2D u grid
276      ! 1/eps
277      unsev(:)=unsev(:)/eps
278      CALL gr_v_scal(1,unsfv,unsfu)
279      CALL gr_v_scal(1,unsev,unseu)
280    END IF
281 
282  END SUBROUTINE ini_slab_transp
283
284  SUBROUTINE divgrad_phy(nlevs,temp,delta)
285! Computes temperature tendency due to horizontal diffusion :
286! T Laplacian, later multiplied by diffusion coef and time-step
287
288    IMPLICIT NONE
289 
290    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
291    REAL, INTENT(IN)   :: temp(klon_glo,nlevs) ! slab temperature
292    REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.)
293    REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs)
294    REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
295    INTEGER :: ll,iip1,jjp1
296 
297    iip1=nbp_lon+1
298    jjp1=nbp_lat
299   
300    ! transpose temp to 2D horiz. grid
301    CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d)
302    ! computes gradient (proportional to heat flx)
303    CALL grad(nlevs,delta_2d,ghx,ghy)
304    ! put flux to 0 at ocean / continent boundary
305    DO ll=1,nlevs
306        ghx(:,ll)=ghx(:,ll)*zmasqu
307        ghy(:,ll)=ghy(:,ll)*zmasqv
308    END DO
309    ! flux divergence
310    CALL diverg(nlevs,ghx,ghy,delta_2d)
311    ! laplacian back to 1D grid
312    CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta)
313 
314    RETURN
315  END SUBROUTINE divgrad_phy
316
317
318  SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy_ek,dt_phy_gm,slab_gm)
319! Temperature tendency for 2-layers slab ocean
320! note : tendency dt later multiplied by (delta time)/(rho.H)
321! to convert from divergence of heat fluxes to T
322
323! 11/16 : Inclusion of GM-like eddy advection
324
325      IMPLICIT NONE
326     
327      LOGICAL,INTENT(in) :: slab_gm
328      ! Here, temperature and flux variables are on 2 layers
329      INTEGER ij
330
331      ! wind stress variables
332      REAL tx_phy(klon_glo),ty_phy(klon_glo)
333      REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
334      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
335      REAL tcurl((nbp_lon+1)*(nbp_lat-1))
336      ! slab temperature on  1D, 2D grid
337      REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
338      ! Temperature gradient, v-points
339      REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat)
340      ! Vertical temperature difference, V-points
341      REAL dtz((nbp_lon+1)*(nbp_lat-1))
342      ! zonal and meridional mass fluxes at u, v points (2D grid)
343      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
344      ! vertical mass flux between the 2 layers
345      REAL fluxv_ek((nbp_lon+1)*nbp_lat)
346      REAL fluxv_gm((nbp_lon+1)*nbp_lat)
347      ! zonal and meridional heat fluxes
348      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
349      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
350      ! temperature tendency (in K.s-1.kg.m-2)
351      REAL dt_ek((nbp_lon+1)*nbp_lat,2), dt_phy_ek(klon_glo,2)
352      REAL dt_gm((nbp_lon+1)*nbp_lat,2), dt_phy_gm(klon_glo,2)
353      ! helper vars
354      REAL zonavg, fluxv
355      REAL, PARAMETER :: sea_den=1025. ! sea water density
356
357      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
358
359! Grid definitions
360      iim=nbp_lon
361      iip1=nbp_lon+1
362      iip2=nbp_lon+2
363      jjp1=nbp_lat
364      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
365      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
366      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
367! Convert temperature to 2D grid
368      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
369
370! ------------------------------------
371! Ekman mass fluxes and Temp tendency
372! ------------------------------------
373! Convert taux,y to 2D  scalar grid
374      ! north and south poles tx,ty no meaning
375      tx_phy(1)=0.
376      tx_phy(klon_glo)=0.
377      ty_phy(1)=0.
378      ty_phy(klon_glo)=0.
379      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
380      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
381
382      IF (ekman_zonavg) THEN ! use zonal average of wind stress
383        DO ij=1,jjp1-2
384          zonavg=SUM(txu(ij*iip1+1:ij*iip1+iim))/iim
385          txu(ij*iip1+1:(ij+1)*iip1)=zonavg
386          zonavg=SUM(tyu(ij*iip1+1:ij*iip1+iim))/iim
387          tyu(ij*iip1+1:(ij+1)*iip1)=zonavg
388        END DO
389      END IF
390
391! ---- Meridional flux ----
392      ! wind stress at v points
393      CALL gr_scal_v(1,txu,txv)
394
395      ! Ekman balance
396      fluxm=-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
397
398      IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
399        tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:) ! dtx/dy
400        !poles curl = 0
401        tcurl(1:iip1)=0.
402        tcurl(ip1jmi1+1:ip1jm)=0.
403        fluxm=fluxm-tcurl*unsbv
404      ENDIF
405
406      IF (slab_tyeq) THEN ! frictional meridional transport at equator
407        CALL gr_scal_v(1,tyu,tyv)
408        fluxm=fluxm+tyv*unsev ! in kg.s-1.m-1 (zonal distance)
409      ENDIF
410
411      !  apply continent mask, multiply by horiz grid dimension
412      fluxm=fluxm*cv*cuvsurcv*zmasqv
413
414! ---- Zonal flux ----
415      IF (ekman_zonadv) THEN
416        CALL gr_scal_u(1,txu,txu) ! wind stress at u points
417        CALL gr_scal_u(1,tyu,tyu)
418        fluxz=tyu*unsfu+txu*unseu ! Ekman + frictional flow
419        !  apply continent mask, multiply by horiz grid dimension
420        fluxz=fluxz*cu*cvusurcu*zmasqu
421      END IF
422           
423!  ---- Vertical mass flux ----
424!  from mass budget (divergence of horiz fluxes)
425      IF (ekman_zonadv) THEN
426         DO ij=iip2,ip1jm
427           fluxv_ek(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
428         ENDDO
429      ELSE
430         DO ij=iip2,ip1jm
431           fluxv_ek(ij)=-fluxm(ij)+fluxm(ij-iip1)
432         ENDDO
433      END IF
434      DO ij=iip1,ip1jmi1,iip1
435         fluxv_ek(ij+1)=fluxv_ek(ij+iip1)
436      END DO
437!  vertical mass flux at Poles
438      fluxv_ek(1)=-SUM(fluxm(1:iim))     
439      fluxv_ek(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
440
441! ---- Meridional heat fluxes ----
442      DO ij=1,ip1jm
443          ! centered scheme
444          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
445          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
446      END DO
447
448! ---- Zonal heat fluxes ----
449! Schema upstream     
450      IF (ekman_zonadv) THEN
451        fluxtz(:,:)=0.     
452        DO ij=iip2,ip1jm
453          IF (fluxz(ij).GE.0.) THEN
454                 fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
455                 fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
456          ELSE
457                 fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
458                 fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
459          ENDIF
460        ENDDO
461        DO ij=iip1,ip1jmp1,iip1
462               fluxtz(ij,:)=fluxtz(ij-iim,:)
463        END DO
464      ELSE
465        fluxtz(:,:)=0.
466      ENDIF       
467                   
468! ---- Temperature tendency, horizontal advection: ----
469      DO ij=iip2,ip1jm
470         dt_ek(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
471                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
472      END DO
473! Poles
474      dt_ek(1,:)=SUM(fluxtm(1:iim,:),dim=1)
475      dt_ek(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
476
477! ------------------------------------
478! GM mass fluxes and Temp tendency
479! ------------------------------------
480     IF (slab_gm) THEN
481! Vertical Temperature difference T1-T2 on v-grid points
482      CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz)
483      dtz(:)=MAX(dtz(:),0.25)
484! Horizontal Temperature differences
485      CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty)
486! Meridional flux = -k.s (s=dyT/dzT)
487! Continent mask, multiply by dz/dy
488      fluxm=dty/dtz*500.*cuvsurcv*zmasqv
489! slope limitation, multiply by kappa
490      fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty)
491! convert to kg/s
492      fluxm(:)=fluxm(:)*sea_den
493! Zonal flux = 0. (temporary)
494      fluxz(:)=0.
495!  Vertical mass flux from mass budget (divergence of horiz fluxes)
496      DO ij=iip2,ip1jm
497        fluxv_gm(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
498      ENDDO
499      DO ij=iip1,ip1jmi1,iip1
500         fluxv_gm(ij+1)=fluxv_gm(ij+iip1)
501      END DO
502!  vertical mass flux at Poles
503      fluxv_gm(1)=-SUM(fluxm(1:iim))     
504      fluxv_gm(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
505
506! Meridional heat fluxes
507      DO ij=1,ip1jm
508          ! centered scheme
509          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
510          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
511      END DO
512
513! Zonal heat fluxes
514! Schema upstream     
515      DO ij=iip2,ip1jm
516      IF (fluxz(ij).GE.0.) THEN
517             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
518             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
519      ELSE
520             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
521             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
522      ENDIF
523      ENDDO
524      DO ij=iip1*2,ip1jmp1,iip1
525             fluxtz(ij,:)=fluxtz(ij-iim,:)
526      END DO
527                   
528! Temperature tendency :
529! divergence of horizontal heat fluxes
530      DO ij=iip2,ip1jm
531         dt_gm(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
532                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
533      END DO
534! Poles
535      dt_gm(1,:)=SUM(fluxtm(1:iim,:),dim=1)
536      dt_gm(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
537     ELSE
538       dt_gm(:,:)=0.
539       fluxv_gm(:)=0.
540     ENDIF ! slab_gm
541
542! ------------------------------------
543! Temp tendency from vertical advection
544! Divide by cell area
545! ------------------------------------
546! vertical heat flux = mass flux * T, upstream scheme
547      DO ij=iip2,ip1jm
548         ! net flux, needed for upstream scheme
549         fluxv=fluxv_ek(ij)+fluxv_gm(ij)
550         ! temperature tendencies from vertical fluxes
551         IF (fluxv.GT.0.) THEN
552           dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,2)
553           dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,2)
554           dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,2)
555           dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,2)
556         ELSE
557           dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,1)
558           dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,1)
559           dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,1)
560           dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,1)
561         ENDIF
562         ! divide by cell area
563         dt_ek(ij,:)=dt_ek(ij,:)/aire(ij)
564         dt_gm(ij,:)=dt_gm(ij,:)/aire(ij)
565      END DO
566      ! North Pole
567      fluxv=fluxv_ek(1)+fluxv_gm(1)
568        IF (fluxv.GT.0.) THEN
569          dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,2)
570          dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,2)
571          dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,2)
572          dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,2)
573        ELSE
574          dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,1)
575          dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,1)
576          dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,1)
577          dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,1)
578        ENDIF
579      dt_ek(1,:)=dt_ek(1,:)/apoln
580      dt_gm(1,:)=dt_gm(1,:)/apoln
581      ! South pole
582      fluxv=fluxv_ek(ip1jmp1)+fluxv_gm(ip1jmp1)
583        IF (fluxv.GT.0.) THEN
584          dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,2)
585          dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,2)
586          dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,2)
587          dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,2)
588        ELSE
589          dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,1)
590          dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,1)
591          dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,1)
592          dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,1)
593        ENDIF
594      dt_ek(ip1jmp1,:)=dt_ek(ip1jmp1,:)/apols
595      dt_gm(ip1jmp1,:)=dt_gm(ip1jmp1,:)/apols
596     
597      dt_ek(2:iip1,1)=dt_ek(1,1)
598      dt_ek(2:iip1,2)=dt_ek(1,2)
599      dt_gm(2:iip1,1)=dt_gm(1,1)
600      dt_gm(2:iip1,2)=dt_gm(1,2)
601      dt_ek(ip1jm+1:ip1jmp1,1)=dt_ek(ip1jmp1,1)
602      dt_ek(ip1jm+1:ip1jmp1,2)=dt_ek(ip1jmp1,2)
603      dt_gm(ip1jm+1:ip1jmp1,1)=dt_gm(ip1jmp1,1)
604      dt_gm(ip1jm+1:ip1jmp1,2)=dt_gm(ip1jmp1,2)
605
606      DO ij=iip1,ip1jmi1,iip1
607         dt_gm(ij+1,:)=dt_gm(ij+iip1,:)
608         dt_ek(ij+1,:)=dt_ek(ij+iip1,:)
609      END DO
610
611! T tendency back to 1D grid...
612      CALL gr_dyn_fi(2,iip1,jjp1,dt_ek,dt_phy_ek)
613      CALL gr_dyn_fi(2,iip1,jjp1,dt_gm,dt_phy_gm)
614
615      RETURN
616  END SUBROUTINE slab_ekman2
617
618  SUBROUTINE slab_gmdiff(ts_phy,dt_phy)
619! Temperature tendency for 2-layers slab ocean
620! Due to Gent-McWilliams type eddy-induced advection
621! Not used anymore ? (included in Ekman)
622     
623      IMPLICIT NONE
624     
625      ! Here, temperature and flux variables are on 2 layers
626      INTEGER ij
627      ! Temperature gradient, v-points
628      REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat)
629      ! Vertical temperature difference, V-points
630      REAL dtz((nbp_lon+1)*(nbp_lat-1))
631      ! slab temperature on  1D, 2D grid
632      REAL ts_phy(klon_glo,2),ts((nbp_lon+1)*nbp_lat,2)
633      ! zonal and meridional mass fluxes at u, v points (2D grid)
634      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
635      ! vertical mass flux between the 2 layers
636      REAL fluxv((nbp_lon+1)*nbp_lat)
637      ! zonal and meridional heat fluxes
638      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
639      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
640      ! temperature tendency (in K.s-1.kg.m-2)
641      REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2)
642
643      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
644
645! Grid definitions
646      iim=nbp_lon
647      iip1=nbp_lon+1
648      iip2=nbp_lon+2
649      jjp1=nbp_lat
650      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
651      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
652      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
653
654! Convert temperature to 2D grid
655      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
656! Vertical Temperature difference T1-T2 on v-grid points
657      CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz)
658      dtz(:)=MAX(dtz(:),0.25)
659! Horizontal Temperature differences
660      CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty)
661! Meridional flux = -k.s (s=dyT/dzT)
662! Continent mask, multiply by dz/dy
663      fluxm=dty/dtz*500.*cuvsurcv*zmasqv
664! slope limitation, multiply by kappa
665      fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty)
666! Zonal flux = 0. (temporary)
667      fluxz(:)=0.
668!  Vertical mass flux from mass budget (divergence of horiz fluxes)
669      DO ij=iip2,ip1jm
670        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
671      ENDDO
672      DO ij=iip1,ip1jmi1,iip1
673         fluxv(ij+1)=fluxv(ij+iip1)
674      END DO
675!  vertical mass flux at Poles
676      fluxv(1)=-SUM(fluxm(1:iim))     
677      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
678      fluxv=fluxv
679
680! Meridional heat fluxes
681      DO ij=1,ip1jm
682          ! centered scheme
683          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
684          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
685      END DO
686
687! Zonal heat fluxes
688! Schema upstream     
689      DO ij=iip2,ip1jm
690      IF (fluxz(ij).GE.0.) THEN
691             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
692             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
693      ELSE
694             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
695             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
696      ENDIF
697      ENDDO
698      DO ij=iip1*2,ip1jmp1,iip1
699             fluxtz(ij,:)=fluxtz(ij-iim,:)
700      END DO
701                   
702! Temperature tendency :
703      DO ij=iip2,ip1jm
704! divergence of horizontal heat fluxes
705         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
706                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
707! + vertical heat flux (mass flux * T, upstream scheme)
708         IF (fluxv(ij).GT.0.) THEN
709           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
710           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
711         ELSE
712           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
713           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
714         ENDIF
715         ! divide by cell area
716         dt(ij,:)=dt(ij,:)/aire(ij)
717      END DO
718      DO ij=iip1,ip1jmi1,iip1
719         dt(ij+1,:)=dt(ij+iip1,:)
720      END DO
721! Poles
722      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
723        IF (fluxv(1).GT.0.) THEN
724          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
725          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
726        ELSE
727          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
728          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
729        ENDIF
730      dt(1,:)=dt(1,:)/apoln
731      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
732       IF (fluxv(ip1jmp1).GT.0.) THEN
733         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
734         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
735       ELSE
736         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
737         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
738       ENDIF
739      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
740      dt(2:iip1,1)=dt(1,1)
741      dt(2:iip1,2)=dt(1,2)
742      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
743      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
744
745! T tendency back to 1D grid...
746      CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy)
747
748      RETURN
749  END SUBROUTINE slab_gmdiff
750
751!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
752!!! Routines copied from dynamics !!!
753!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
754
755  SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn)
756  ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
757  IMPLICIT NONE
758
759  INTEGER,INTENT(IN) :: im,jm,nfield
760  REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid
761  REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid
762
763  INTEGER :: i,j,ifield,ig
764
765  DO ifield=1,nfield
766    ! Handle poles
767    DO i=1,im
768      pdyn(i,1,ifield)=pfi(1,ifield)
769      pdyn(i,jm,ifield)=pfi(klon_glo,ifield)
770    ENDDO
771    ! Other points
772    DO j=2,jm-1
773      ig=2+(j-2)*(im-1)
774      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
775      pdyn(im,j,ifield)=pdyn(1,j,ifield)
776    ENDDO
777  ENDDO ! of DO ifield=1,nfield
778
779  END SUBROUTINE gr_fi_dyn
780
781  SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi)
782  ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
783  IMPLICIT NONE
784
785  INTEGER,INTENT(IN) :: im,jm,nfield
786  REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid
787  REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid
788
789  INTEGER j,ifield,ig
790
791  ! Sanity check:
792  IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN
793    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
794    STOP
795  ENDIF
796
797  ! Handle poles
798  CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo)
799  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo)
800  ! Other points
801  DO ifield=1,nfield
802    DO j=2,jm-1
803      ig=2+(j-2)*(im-1)
804      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
805    ENDDO
806  ENDDO
807
808  END SUBROUTINE gr_dyn_fi
809
810  SUBROUTINE  grad(klevel,pg,pgx,pgy)
811  ! compute the covariant components pgx,pgy of the gradient of pg
812  ! pgx = d(pg)/dx * delta(x) = delta(pg)
813  IMPLICIT NONE
814
815  INTEGER,INTENT(IN) :: klevel
816  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
817  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
818  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
819
820  INTEGER :: l,ij
821  INTEGER :: iim,iip1,ip1jm,ip1jmp1
822
823  iim=nbp_lon
824  iip1=nbp_lon+1
825  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
826  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
827
828  DO l=1,klevel
829    DO ij=1,ip1jmp1-1
830      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
831    ENDDO
832    ! correction for pgx(ip1,j,l) ...
833    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
834    DO ij=iip1,ip1jmp1,iip1
835      pgx(ij,l)=pgx(ij-iim,l)
836    ENDDO
837    DO ij=1,ip1jm
838      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
839    ENDDO
840  ENDDO
841
842  END SUBROUTINE grad
843
844  SUBROUTINE diverg(klevel,x,y,div)
845  ! computes the divergence of a vector field of components
846  ! x,y. x and y being covariant components
847  IMPLICIT NONE
848
849  INTEGER,INTENT(IN) :: klevel
850  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
851  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
852  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
853
854  INTEGER :: l,ij
855  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
856
857  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
858  REAL :: sumypn,sumyps
859  REAL,EXTERNAL :: SSUM
860
861  iim=nbp_lon
862  iip1=nbp_lon+1
863  iip2=nbp_lon+2
864  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
865  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
866  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
867
868  DO l=1,klevel
869    DO ij=iip2,ip1jm-1
870      div(ij+1,l)= &
871        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
872        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
873    ENDDO
874    ! correction for div(1,j,l) ...
875    ! ... div(1,j,l)= div(iip1,j,l) ...
876    DO ij=iip2,ip1jm,iip1
877      div(ij,l)=div(ij+iim,l)
878    ENDDO
879    ! at the poles
880    DO ij=1,iim
881      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
882      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
883    ENDDO
884    sumypn=SSUM(iim,aiy1,1)/apoln
885    sumyps=SSUM(iim,aiy2,1)/apols
886    DO ij=1,iip1
887      div(ij,l)=-sumypn
888      div(ij+ip1jm,l)=sumyps
889    ENDDO
890    ! End (poles)
891  ENDDO ! of DO l=1,klevel
892
893  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
894  DO l=1,klevel
895    DO ij=iip2,ip1jm
896      div(ij,l)=div(ij,l)*unsaire(ij)
897    ENDDO
898  ENDDO
899
900  END SUBROUTINE diverg
901
902  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
903  ! convert values from v points to scalar points on C-grid
904  ! used to  compute unsfu, unseu (u points, but depends only on latitude)
905  IMPLICIT NONE
906
907  INTEGER,INTENT(IN) :: nx ! number of levels or fields
908  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
909  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
910
911  INTEGER :: l,ij
912  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
913
914  iip1=nbp_lon+1
915  iip2=nbp_lon+2
916  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
917  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
918
919  DO l=1,nx
920    DO ij=iip2,ip1jm
921      x_scal(ij,l)= &
922                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
923                  /(airev(ij-iip1)+airev(ij))
924    ENDDO
925    DO ij=1,iip1
926      x_scal(ij,l)=0.
927    ENDDO
928    DO ij=ip1jm+1,ip1jmp1
929      x_scal(ij,l)=0.
930    ENDDO
931  ENDDO
932
933  END SUBROUTINE gr_v_scal
934 
935  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
936  ! convert values from scalar points to v points on C-grid
937  ! used to compute wind stress at V points
938  IMPLICIT NONE
939
940  INTEGER,INTENT(IN) :: nx ! number of levels or fields
941  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
942  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
943
944  INTEGER :: l,ij
945  INTEGER :: iip1,ip1jm
946
947  iip1=nbp_lon+1
948  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
949
950      DO l=1,nx
951        DO ij=1,ip1jm
952          x_v(ij,l)= &
953            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
954            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
955            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
956        ENDDO
957      ENDDO
958
959  END SUBROUTINE gr_scal_v
960
961  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
962  ! convert values from scalar points to U points on C-grid
963  ! used to compute wind stress at U points
964  IMPLICIT NONE
965
966  INTEGER,INTENT(IN) :: nx
967  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
968  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
969
970  INTEGER :: l,ij
971  INTEGER :: iip1,jjp1,ip1jmp1
972
973  iip1=nbp_lon+1
974  jjp1=nbp_lat
975  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
976
977  DO l=1,nx
978     DO ij=1,ip1jmp1-1
979        x_u(ij,l)= &
980         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
981         /(aire(ij)+aire(ij+1))
982     ENDDO
983  ENDDO
984
985  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
986
987  END SUBROUTINE gr_scal_u
988
989END MODULE slab_heat_transp_mod
Note: See TracBrowser for help on using the repository browser.