source: LMDZ6/trunk/libf/phylmd/slab_heat_transp_mod.F90 @ 3504

Last change on this file since 3504 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

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