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

Last change on this file since 4350 was 3531, checked in by Laurent Fairhead, 5 years ago

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