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

Last change on this file since 2793 was 2656, checked in by Ehouarn Millour, 8 years ago

Making the slab work:

  • added a slab_heat_transp_mod module for horizontal diffusion and Ekman transport
  • added storage and output of relevent variables in phyredem, phyetat0, phy_output_ctrlout_mod, phys_output_write_mod
  • moved nslay (number of slab layers) out of dimphy into ocean_slab_mod.

FC
File size: 23.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 :: unsairez(:) ! 1/cell area
14  !$OMP THREADPRIVATE(unsairez)
15  REAL,SAVE,ALLOCATABLE :: unsaire(:)
16  !$OMP THREADPRIVATE(unsaire)
17  REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m)
18  !$OMP THREADPRIVATE(cu)
19  REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m)
20  !$OMP THREADPRIVATE(cv)
21  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points)
22  !$OMP THREADPRIVATE(cuvsurcv)
23  REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points)
24  !$OMP THREADPRIVATE(cvusurcu)
25  REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area
26  !$OMP THREADPRIVATE(aire)
27  REAL,SAVE :: apoln ! area of north pole points
28  !$OMP THREADPRIVATE(apoln)
29  REAL,SAVE :: apols ! area of south pole points
30  !$OMP THREADPRIVATE(apols)
31  REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells
32  !$OMP THREADPRIVATE(aireu)
33  REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells
34  !$OMP THREADPRIVATE(airev)
35
36  ! Local variables for horiz mass flux in slab
37  LOGICAL,SAVE :: alpha_var
38  !$OMP THREADPRIVATE(alpha_var)
39  LOGICAL,SAVE :: slab_upstream
40  !$OMP THREADPRIVATE(slab_upstream)
41
42  REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux
43  !$OMP THREADPRIVATE(zmasqu)
44  REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux
45  !$OMP THREADPRIVATE(zmasqv)
46  REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points
47  !$OMP THREADPRIVATE(unsfv)
48  REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag)
49  !$OMP THREADPRIVATE(unsev)
50  REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points
51  !$OMP THREADPRIVATE(unsfu)
52  REAL,SAVE,ALLOCATABLE :: unseu(:)
53  !$OMP THREADPRIVATE(unseu)
54
55  ! Routines from dyn3d, valid on global dynamics grid only:
56  PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid
57  PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points
58  PRIVATE :: grad,diverg
59
60CONTAINS
61 
62  SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
63                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
64                                  aire_,apoln_,apols_, &
65                                  aireu_,airev_)
66    ! number of points in lon, lat
67    IMPLICIT NONE
68    ! Routine copies some geometry variables from the dynamical core
69    ! see global vars for meaning
70    INTEGER,INTENT(IN) :: ip1jm
71    INTEGER,INTENT(IN) :: ip1jmp1
72    REAL,INTENT(IN) :: unsairez_(ip1jm)
73    REAL,INTENT(IN) :: fext_(ip1jm)
74    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
75    REAL,INTENT(IN) :: cu_(ip1jmp1)
76    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
77    REAL,INTENT(IN) :: cv_(ip1jm)
78    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
79    REAL,INTENT(IN) :: aire_(ip1jmp1)
80    REAL,INTENT(IN) :: apoln_
81    REAL,INTENT(IN) :: apols_
82    REAL,INTENT(IN) :: aireu_(ip1jmp1)
83    REAL,INTENT(IN) :: airev_(ip1jm)
84
85
86    ! Sanity check on dimensions
87    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
88        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
89      write(*,*) "ini_slab_transp_geom Error: wrong array sizes"
90      stop
91    endif
92! Allocations could be done only on master process/thread...
93    allocate(unsairez(ip1jm))
94    unsairez(:)=unsairez_(:)
95    allocate(fext(ip1jm))
96    fext(:)=fext_(:)
97    allocate(unsaire(ip1jmp1))
98    unsaire(:)=unsaire_(:)
99    allocate(cu(ip1jmp1))
100    cu(:)=cu_(:)
101    allocate(cuvsurcv(ip1jm))
102    cuvsurcv(:)=cuvsurcv_(:)
103    allocate(cv(ip1jm))
104    cv(:)=cv_(:)
105    allocate(cvusurcu(ip1jmp1))
106    cvusurcu(:)=cvusurcu_(:)
107    allocate(aire(ip1jmp1))
108    aire(:)=aire_(:)
109    apoln=apoln_
110    apols=apols_
111    allocate(aireu(ip1jmp1))
112    aireu(:)=aireu_(:)
113    allocate(airev(ip1jm))
114    airev(:)=airev_(:)
115
116  END SUBROUTINE ini_slab_transp_geom
117
118  SUBROUTINE ini_slab_transp(zmasq)
119
120!    USE ioipsl_getin_p_mod, only: getin_p
121    USE IOIPSL, ONLY : getin
122    IMPLICIT NONE
123
124    REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent
125    REAL zmasq_2d((nbp_lon+1)*nbp_lat)
126    REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter
127    REAL eps ! epsilon friction timescale (s-1)
128    INTEGER :: slab_ekman
129    INTEGER i
130    INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
131
132! Some definition for grid size
133    ip1jm=(nbp_lon+1)*(nbp_lat-1)
134    ip1jmp1=(nbp_lon+1)*nbp_lat
135    iim=nbp_lon
136    iip1=nbp_lon+1
137    jjp1=nbp_lat
138    ip1jm=(nbp_lon+1)*(nbp_lat-1)
139    ip1jmp1=(nbp_lon+1)*nbp_lat
140
141! Define ocean / continent mask (no flux into continent cell)
142    allocate(zmasqu(ip1jmp1))
143    allocate(zmasqv(ip1jm))
144    zmasqu=1.
145    zmasqv=1.
146
147    ! convert mask to 2D grid
148    CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d)
149    ! put flux mask to 0 at boundaries of continent cells
150    DO i=1,ip1jmp1-1
151      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
152              zmasqu(i)=0.
153      ENDIF
154    END DO
155    DO i=iip1,ip1jmp1,iip1
156      zmasqu(i)=zmasqu(i-iim)
157    END DO
158    DO i=1,ip1jm
159      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN
160              zmasqv(i)=0.
161      END IF
162    END DO
163    slab_ekman=2
164    CALL getin("slab_ekman",slab_ekman)
165! Coriolis and friction for Ekman transport
166    IF (slab_ekman.GT.0) THEN
167      allocate(unsfv(ip1jm))
168      allocate(unsev(ip1jm))
169      allocate(unsfu(ip1jmp1))
170      allocate(unseu(ip1jmp1))
171
172      eps=1e-5 ! Drag
173      CALL getin('slab_eps',eps)
174      print *,'epsilon=',eps
175      ff=fext*unsairez ! Coriolis
176      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
177      ! on 2D grid v points
178      DO i=1,ip1jm
179        unsev(i)=eps/(ff(i)*ff(i)+eps**2)
180        unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2)
181      ENDDO
182      ! compute values on 2D u grid
183      CALL gr_v_scal(1,unsfv,unsfu)
184      CALL gr_v_scal(1,unsev,unseu)
185    END IF
186
187! Other options for Ekman transport
188    ! Alpha variable?
189      alpha_var=.FALSE.
190      CALL getin('slab_alphav',alpha_var)
191      print *,'alpha variable',alpha_var
192!  centered ou upstream scheme for meridional transport
193      slab_upstream=.FALSE.
194      CALL getin('slab_upstream',slab_upstream)
195      print *,'upstream slab scheme',slab_upstream
196 
197  END SUBROUTINE ini_slab_transp
198
199  SUBROUTINE divgrad_phy(nlevs,temp,delta)
200! Computes temperature tendency due to horizontal diffusion :
201! T Laplacian, later multiplied by diffusion coef and time-step
202
203    IMPLICIT NONE
204 
205    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
206    REAL, INTENT(IN)   :: temp(klon_glo,nlevs) ! slab temperature
207    REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.)
208    REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs)
209    REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
210    INTEGER :: ll,iip1,jjp1
211 
212    iip1=nbp_lon+1
213    jjp1=nbp_lat
214   
215    ! transpose temp to 2D horiz. grid
216    CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d)
217    ! computes gradient (proportional to heat flx)
218    CALL grad(nlevs,delta_2d,ghx,ghy)
219    ! put flux to 0 at ocean / continent boundary
220    DO ll=1,nlevs
221        ghx(:,ll)=ghx(:,ll)*zmasqu
222        ghy(:,ll)=ghy(:,ll)*zmasqv
223    END DO
224    ! flux divergence
225    CALL diverg(nlevs,ghx,ghy,delta_2d)
226    ! laplacian back to 1D grid
227    CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta)
228 
229    RETURN
230  END SUBROUTINE divgrad_phy
231
232  SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy)
233! 1.5 Layer Ekman transport temperature tendency
234! note : tendency dt later multiplied by (delta t)/(rho.H)
235! to convert from divergence of heat fluxes to T
236
237      IMPLICIT NONE
238     
239      ! tx, ty : wind stress (different grids)
240      ! fluxm, fluz : mass *or heat* fluxes
241      ! dt : temperature tendency
242      INTEGER ij
243
244      ! ts surface temp, td deep temp (diagnosed)
245      REAL ts_phy(klon_glo)
246      REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat)
247      ! zonal and meridional wind stress
248      REAL tx_phy(klon_glo),ty_phy(klon_glo)
249      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
250      REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1))
251      ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
252      REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1))
253      ! vertical  and absolute mass fluxes (to estimate alpha)
254      REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1))
255      ! temperature tendency
256      REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo)
257      REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef
258
259      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
260
261! Grid definitions
262      iim=nbp_lon
263      iip1=nbp_lon+1
264      iip2=nbp_lon+2
265      jjp1=nbp_lat
266      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
267      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
268      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
269
270! Convert taux,y to 2D  scalar grid
271! Note: 2D grid size = iim*jjm. iip1=iim+1
272! First and last points in zonal direction are the same
273! we use 1 index ij from 1 to (iim+1)*(jjm+1)
274      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
275      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
276! convert to u,v grid (Arakawa C)
277! Multiply by f or eps to get mass flux
278      ! Meridional mass flux
279      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
280      CALL gr_scal_v(1,tyu,tyv)
281      fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
282      ! Zonal mass flux
283      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
284      CALL gr_scal_u(1,tyu,tyu)
285      fluxz=tyu*unsfu+txu*unseu
286           
287! Correct flux: continent mask and horiz grid size
288      ! multiply m-flux by mask and dx: flux in kg.s-1
289      fluxm=fluxm*cv*cuvsurcv*zmasqv
290      ! multiply z-flux by mask and dy: flux in kg.s-1
291      fluxz=fluxz*cu*cvusurcu*zmasqu
292
293! Compute vertical  and absolute mass flux (for variable alpha)
294      IF (alpha_var) THEN
295        DO ij=iip2,ip1jm
296        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
297        fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) &
298               +ABS(fluxm(ij))+ABS(fluxm(ij-iip1))
299        ENDDO
300        DO ij=iip1,ip1jmi1,iip1
301            fluxt(ij+1)=fluxt(ij+iip1)
302            fluxv(ij+1)=fluxv(ij+iip1)
303        END DO
304        fluxt(1)=SUM(ABS(fluxm(1:iim)))
305        fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1)))
306        fluxv(1)=-SUM(fluxm(1:iim))
307        fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
308        fluxt=MAX(fluxt,1.e-10)
309      ENDIF
310
311! Compute alpha coefficient.
312! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
313      IF (alpha_var) THEN
314          ! increase alpha (and Tdeep) in downwelling regions
315          ! and decrease in upwelling regions
316          ! to avoid "hot spots" where there is surface convergence
317          DO ij=iip2,ip1jm
318              alpha(ij)=2./3.-fluxv(ij)/fluxt(ij)/3.
319          ENDDO
320          alpha(1:iip1)=2./3.-fluxv(1)/fluxt(1)/3.
321          alpha(ip1jm+1:ip1jmp1)=2./3.-fluxv(ip1jmp1)/fluxt(ip1jmp1)/3.
322      ELSE
323          alpha(:)=2./3.
324          ! Tsurf-Tdeep ~ 10° in the Tropics
325      ENDIF
326
327! Estimate deep temperature
328      CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts)
329      DO ij=1,ip1jmp1
330         td(ij)=271.35+(ts(ij)-271.35)*alpha(ij)
331         td(ij)=MIN(td(ij),ts(ij))
332      END DO
333       
334! Meridional heat flux: multiply mass flux by (ts-td)
335! flux in K.kg.s-1
336      IF (slab_upstream) THEN
337        ! upstream scheme to avoid hot spots
338        DO ij=1,ip1jm
339        IF (fluxm(ij).GE.0.) THEN
340           fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij))
341        ELSE
342           fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1))
343        END IF
344        END DO
345      ELSE
346        ! centered scheme better in mid-latitudes
347        DO ij=1,ip1jm
348          fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2.
349        END DO
350      ENDIF
351
352! Zonal heat flux
353! upstream scheme
354      DO ij=iip2,ip1jm
355          fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2.
356      END DO
357      DO ij=iip1*2,ip1jmp1,iip1
358           fluxz(ij)=fluxz(ij-iim)
359      END DO
360
361! temperature tendency = divergence of heat fluxes
362! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
363      DO ij=iip2,ip1jm
364         dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) &
365               /aire(ij) ! aire : grid area
366      END DO
367      DO ij=iip1,ip1jmi1,iip1
368         dt(ij+1)=dt(ij+iip1)
369      END DO
370! special treatment at the Poles
371      dt(1)=SUM(fluxm(1:iim))/apoln
372      dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols
373      dt(2:iip1)=dt(1)
374      dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1)
375
376! tendencies back to 1D grid
377      CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy)
378
379      RETURN
380  END SUBROUTINE slab_ekman1
381
382  SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy)
383! Temperature tendency for 2-layers slab ocean
384! note : tendency dt later multiplied by (delta time)/(rho.H)
385! to convert from divergence of heat fluxes to T
386     
387      IMPLICIT NONE
388     
389      ! Here, temperature and flux variables are on 2 layers
390      INTEGER ij
391
392      REAL tx_phy(klon_glo),ty_phy(klon_glo)
393      REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
394      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
395      ! slab temperature on  1D, 2D grid
396      REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
397      ! zonal and meridional Ekman mass flux at u, v points (2D grid)
398      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
399      ! vertical mass flux between the 2 layers
400      REAL fluxv((nbp_lon+1)*nbp_lat)
401      ! zonal and meridional heat fluxes
402      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
403      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
404      ! temperature tendency (in K.s-1.kg.m-2)
405      REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2)
406
407      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
408
409! Grid definitions
410      iim=nbp_lon
411      iip1=nbp_lon+1
412      iip2=nbp_lon+2
413      jjp1=nbp_lat
414      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
415      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
416      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
417
418! Convert taux,y to 2D  scalar grid
419      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
420      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
421! Multiply taux,y by f or eps, and convert to 2D u,v grids
422! (Arakawa C grid)
423      ! Meridional flux
424      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
425      CALL gr_scal_v(1,tyu,tyv)
426      fluxm=tyv*unsev-txv*unsfv
427      ! Zonal flux
428      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
429      CALL gr_scal_u(1,tyu,tyu)
430      fluxz=tyu*unsfu+txu*unseu
431           
432! Convert temperature to 2D grid
433      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
434       
435!  Mass fluxes (apply continent mask, multiply by horiz grid dimension)
436      fluxm=fluxm*cv*cuvsurcv*zmasqv
437      fluxz=fluxz*cu*cvusurcu*zmasqu
438!  Vertical mass flux from mass budget (divergence of horiz fluxes)
439      DO ij=iip2,ip1jm
440        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
441      ENDDO
442      DO ij=iip1,ip1jmi1,iip1
443         fluxv(ij+1)=fluxv(ij+iip1)
444      END DO
445!  vertical mass flux at Poles
446      fluxv(1)=-SUM(fluxm(1:iim))     
447      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
448      fluxv=fluxv
449
450! Meridional heat fluxes
451      DO ij=1,ip1jm
452          ! centered scheme
453          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
454          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
455      END DO
456
457! Zonal heat fluxes
458! Schema upstream     
459      DO ij=iip2,ip1jm
460      IF (fluxz(ij).GE.0.) THEN
461             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
462             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
463      ELSE
464             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
465             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
466      ENDIF
467      ENDDO
468      DO ij=iip1*2,ip1jmp1,iip1
469             fluxtz(ij,:)=fluxtz(ij-iim,:)
470      END DO
471                   
472! Temperature tendency :
473      DO ij=iip2,ip1jm
474! divergence of horizontal heat fluxes
475         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
476                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
477! + vertical heat flux (mass flux * T, upstream scheme)
478         IF (fluxv(ij).GT.0.) THEN
479           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
480           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
481         ELSE
482           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
483           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
484         ENDIF
485         ! divide by cell area
486         dt(ij,:)=dt(ij,:)/aire(ij)
487      END DO
488      DO ij=iip1,ip1jmi1,iip1
489         dt(ij+1,:)=dt(ij+iip1,:)
490      END DO
491! Pôles
492      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
493        IF (fluxv(1).GT.0.) THEN
494          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
495          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
496        ELSE
497          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
498          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
499        ENDIF
500      dt(1,:)=dt(1,:)/apoln
501      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
502       IF (fluxv(ip1jmp1).GT.0.) THEN
503         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
504         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
505       ELSE
506         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
507         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
508       ENDIF
509      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
510      dt(2:iip1,1)=dt(1,1)
511      dt(2:iip1,2)=dt(1,2)
512      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
513      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
514
515! T tendency back to 1D grid...
516      CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy)
517
518      RETURN
519  END SUBROUTINE slab_ekman2
520
521!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522
523  SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn)
524  ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
525  IMPLICIT NONE
526
527  INTEGER,INTENT(IN) :: im,jm,nfield
528  REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid
529  REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid
530
531  INTEGER :: i,j,ifield,ig
532
533  DO ifield=1,nfield
534    ! Handle poles
535    DO i=1,im
536      pdyn(i,1,ifield)=pfi(1,ifield)
537      pdyn(i,jm,ifield)=pfi(klon_glo,ifield)
538    ENDDO
539    ! Other points
540    DO j=2,jm-1
541      ig=2+(j-2)*(im-1)
542      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
543      pdyn(im,j,ifield)=pdyn(1,j,ifield)
544    ENDDO
545  ENDDO ! of DO ifield=1,nfield
546
547  END SUBROUTINE gr_fi_dyn
548
549  SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi)
550  ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
551  IMPLICIT NONE
552
553  INTEGER,INTENT(IN) :: im,jm,nfield
554  REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid
555  REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid
556
557  INTEGER j,ifield,ig
558
559  ! Sanity check:
560  IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN
561    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
562    STOP
563  ENDIF
564
565  ! Handle poles
566  CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo)
567  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo)
568  ! Other points
569  DO ifield=1,nfield
570    DO j=2,jm-1
571      ig=2+(j-2)*(im-1)
572      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
573    ENDDO
574  ENDDO
575
576  END SUBROUTINE gr_dyn_fi
577
578  SUBROUTINE  grad(klevel,pg,pgx,pgy)
579  ! compute the covariant components pgx,pgy of the gradient of pg
580  ! pgx = d(pg)/dx * delta(x) = delta(pg)
581  IMPLICIT NONE
582
583  INTEGER,INTENT(IN) :: klevel
584  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
585  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
586  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
587
588  INTEGER :: l,ij
589  INTEGER :: iim,iip1,ip1jm,ip1jmp1
590
591  iim=nbp_lon
592  iip1=nbp_lon+1
593  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
594  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
595
596  DO l=1,klevel
597    DO ij=1,ip1jmp1-1
598      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
599    ENDDO
600    ! correction for pgx(ip1,j,l) ...
601    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
602    DO ij=iip1,ip1jmp1,iip1
603      pgx(ij,l)=pgx(ij-iim,l)
604    ENDDO
605    DO ij=1,ip1jm
606      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
607    ENDDO
608  ENDDO
609
610  END SUBROUTINE grad
611
612  SUBROUTINE diverg(klevel,x,y,div)
613  ! computes the divergence of a vector field of components
614  ! x,y. x and y being covariant components
615  IMPLICIT NONE
616
617  INTEGER,INTENT(IN) :: klevel
618  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
619  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
620  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
621
622  INTEGER :: l,ij
623  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
624
625  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
626  REAL :: sumypn,sumyps
627  REAL,EXTERNAL :: SSUM
628
629  iim=nbp_lon
630  iip1=nbp_lon+1
631  iip2=nbp_lon+2
632  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
633  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
634  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
635
636  DO l=1,klevel
637    DO ij=iip2,ip1jm-1
638      div(ij+1,l)= &
639        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
640        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
641    ENDDO
642    ! correction for div(1,j,l) ...
643    ! ... div(1,j,l)= div(iip1,j,l) ...
644    DO ij=iip2,ip1jm,iip1
645      div(ij,l)=div(ij+iim,l)
646    ENDDO
647    ! at the poles
648    DO ij=1,iim
649      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
650      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
651    ENDDO
652    sumypn=SSUM(iim,aiy1,1)/apoln
653    sumyps=SSUM(iim,aiy2,1)/apols
654    DO ij=1,iip1
655      div(ij,l)=-sumypn
656      div(ij+ip1jm,l)=sumyps
657    ENDDO
658    ! End (poles)
659  ENDDO ! of DO l=1,klevel
660
661  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
662  DO l=1,klevel
663    DO ij=iip2,ip1jm
664      div(ij,l)=div(ij,l)*unsaire(ij)
665    ENDDO
666  ENDDO
667
668  END SUBROUTINE diverg
669
670  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
671  ! convert values from v points to scalar points on C-grid
672  ! used to  compute unsfu, unseu (u points, but depends only on latitude)
673  IMPLICIT NONE
674
675  INTEGER,INTENT(IN) :: nx ! number of levels or fields
676  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
677  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
678
679  INTEGER :: l,ij
680  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
681
682  iip1=nbp_lon+1
683  iip2=nbp_lon+2
684  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
685  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
686
687  DO l=1,nx
688    DO ij=iip2,ip1jm
689      x_scal(ij,l)= &
690                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
691                  /(airev(ij-iip1)+airev(ij))
692    ENDDO
693    DO ij=1,iip1
694      x_scal(ij,l)=0.
695    ENDDO
696    DO ij=ip1jm+1,ip1jmp1
697      x_scal(ij,l)=0.
698    ENDDO
699  ENDDO
700
701  END SUBROUTINE gr_v_scal
702 
703  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
704  ! convert values from scalar points to v points on C-grid
705  ! used to compute wind stress at V points
706  IMPLICIT NONE
707
708  INTEGER,INTENT(IN) :: nx ! number of levels or fields
709  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
710  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
711
712  INTEGER :: l,ij
713  INTEGER :: iip1,ip1jm
714
715  iip1=nbp_lon+1
716  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
717
718      DO l=1,nx
719        DO ij=1,ip1jm
720          x_v(ij,l)= &
721            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
722            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
723            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
724        ENDDO
725      ENDDO
726
727  END SUBROUTINE gr_scal_v
728
729  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
730  ! convert values from scalar points to U points on C-grid
731  ! used to compute wind stress at U points
732  IMPLICIT NONE
733
734  INTEGER,INTENT(IN) :: nx
735  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
736  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
737
738  INTEGER :: l,ij
739  INTEGER :: iip1,jjp1,ip1jmp1
740
741  iip1=nbp_lon+1
742  jjp1=nbp_lat
743  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
744
745  DO l=1,nx
746     DO ij=1,ip1jmp1-1
747        x_u(ij,l)= &
748         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
749         /(aire(ij)+aire(ij+1))
750     ENDDO
751  ENDDO
752
753  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
754
755  END SUBROUTINE gr_scal_u
756
757END MODULE slab_heat_transp_mod
Note: See TracBrowser for help on using the repository browser.