source: trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90 @ 2987

Last change on this file since 2987 was 2111, checked in by emillour, 6 years ago

Generic model:

  • some fixes for the slab ocean. Still need to make it work in parallel.

EM
  • Property svn:executable set to *
File size: 15.7 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ocean_slab_mod.F90,v 1.3 2008-02-04 16:24:28 fairhead Exp $
3!
4MODULE surf_heat_transp_mod
5
6IMPLICIT NONE 
7
8  ! Variables copied over from dyn3d dynamics:
9  REAL,SAVE,ALLOCATABLE :: fext(:)
10  REAL,SAVE,ALLOCATABLE :: unsairez(:)
11  REAL,SAVE,ALLOCATABLE :: unsaire(:)
12  REAL,SAVE,ALLOCATABLE :: cu(:)
13  REAL,SAVE,ALLOCATABLE :: cv(:)
14  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:)
15  REAL,SAVE,ALLOCATABLE :: cvusurcu(:)
16  REAL,SAVE,ALLOCATABLE :: aire(:)
17  REAL,SAVE :: apoln
18  REAL,SAVE :: apols
19  REAL,SAVE,ALLOCATABLE :: aireu(:)
20  REAL,SAVE,ALLOCATABLE :: airev(:)
21
22  LOGICAL,SAVE :: alpha_var
23  LOGICAL,SAVE :: slab_upstream
24  REAL,SAVE,ALLOCATABLE :: zmasqu(:)
25  REAL,SAVE,ALLOCATABLE :: zmasqv(:)
26  REAL,SAVE,ALLOCATABLE :: unsfv(:)
27  REAL,SAVE,ALLOCATABLE :: unsev(:)
28  REAL,SAVE,ALLOCATABLE :: unsfu(:)
29  REAL,SAVE,ALLOCATABLE :: unseu(:)
30
31  ! Routines usable only by routines within this module:
32  PRIVATE :: gr_fi_dyn, gr_dyn_fi
33  ! Routines from dyn3d, valid on global dynamics grid only:
34  PRIVATE :: grad,diverg,gr_v_scal,gr_scal_v,gr_scal_u
35 
36CONTAINS
37 
38  SUBROUTINE ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
39                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
40                                  aire_,apoln_,apols_, &
41                                  aireu_,airev_)
42    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
43    IMPLICIT NONE
44    ! Transfer some variables from dyn3d dynamics
45    INTEGER,INTENT(IN) :: ip1jm
46    INTEGER,INTENT(IN) :: ip1jmp1
47    REAL,INTENT(IN) :: unsairez_(ip1jm)
48    REAL,INTENT(IN) :: fext_(ip1jm)
49    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
50    REAL,INTENT(IN) :: cu_(ip1jmp1)
51    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
52    REAL,INTENT(IN) :: cv_(ip1jm)
53    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
54    REAL,INTENT(IN) :: aire_(ip1jmp1)
55    REAL,INTENT(IN) :: apoln_
56    REAL,INTENT(IN) :: apols_
57    REAL,INTENT(IN) :: aireu_(ip1jmp1)
58    REAL,INTENT(IN) :: airev_(ip1jm)
59   
60    ! Sanity check
61    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
62        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
63      write(*,*) "ini_surf_heat_transp Error: wrong array sizes"
64      stop
65    endif
66   
67    allocate(unsairez(ip1jm))
68    unsairez(:)=unsairez_(:)
69    allocate(fext(ip1jm))
70    fext(:)=fext_(:)
71    allocate(unsaire(ip1jmp1))
72    unsaire(:)=unsaire_(:)
73    allocate(cu(ip1jmp1))
74    cu(:)=cu_(:)
75    allocate(cuvsurcv(ip1jm))
76    cuvsurcv(:)=cuvsurcv_(:)
77    allocate(cv(ip1jm))
78    cv(:)=cv_(:)
79    allocate(cvusurcu(ip1jmp1))
80    cvusurcu(:)=cvusurcu_(:)
81    allocate(aire(ip1jmp1))
82    aire(:)=aire_(:)
83    apoln=apoln_
84    apols=apols_
85    allocate(aireu(ip1jmp1))
86    aireu(:)=aireu_(:)
87    allocate(airev(ip1jm))
88    airev(:)=airev_(:)
89   
90  END SUBROUTINE ini_surf_heat_transp
91 
92  SUBROUTINE ini_surf_heat_transp_mod
93    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
94    IMPLICIT NONE
95    INTEGER :: ip1jm, ip1jmp1
96   
97    ip1jm=(nbp_lon+1)*(nbp_lat-1)
98    ip1jmp1=(nbp_lon+1)*nbp_lat
99 
100    allocate(zmasqu(ip1jmp1))
101    allocate(zmasqv(ip1jm))
102    allocate(unsfv(ip1jm))
103    allocate(unsev(ip1jm))
104    allocate(unsfu(ip1jmp1))
105    allocate(unseu(ip1jmp1))
106
107  END SUBROUTINE ini_surf_heat_transp_mod
108
109  SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
110
111      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
112      IMPLICIT NONE
113
114      INTEGER,INTENT(IN) :: ngrid, nlevs
115      REAL,INTENT(IN) :: temp(ngrid,nlevs)
116      REAL,INTENT(OUT) :: delta(ngrid,nlevs)
117      REAL delta_2d((nbp_lon+1)*nbp_lat,nlevs)
118      INTEGER :: ll
119      REAL ghx((nbp_lon+1)*nbp_lat,nlevs)
120      REAL ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
121      INTEGER :: iip1,jjp1
122     
123      iip1=nbp_lon+1
124      jjp1=nbp_lat
125
126      CALL gr_fi_dyn(nlevs,ngrid,iip1,jjp1,temp,delta_2d)
127      CALL grad(nlevs,delta_2d,ghx,ghy)
128      DO ll=1,nlevs
129          ghx(:,ll)=ghx(:,ll)*zmasqu
130! pas de diffusion zonale 
131!          ghx(:,ll)=0.
132          ghy(:,ll)=ghy(:,ll)*zmasqv
133      END DO
134      CALL diverg(nlevs,ghx,ghy,delta_2d)
135      CALL gr_dyn_fi(nlevs,iip1,jjp1,ngrid,delta_2d,delta)
136
137
138  END SUBROUTINE divgrad_phy
139
140
141
142  SUBROUTINE init_masquv(ngrid,zmasq)
143
144      USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
145      IMPLICIT NONE
146
147
148      INTEGER,INTENT(IN) :: ngrid
149      REAL zmasq(ngrid)
150      REAL zmasq_2d((nbp_lon+1)*nbp_lat)
151      REAL ff((nbp_lon+1)*(nbp_lat-1))
152      REAL eps
153      INTEGER i
154      INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
155     
156      iim=nbp_lon
157      iip1=nbp_lon+1
158      jjp1=nbp_lat
159      ip1jm=(nbp_lon+1)*(nbp_lat-1)
160      ip1jmp1=(nbp_lon+1)*nbp_lat
161
162! Masques u,v
163      zmasqu=1.
164      zmasqv=1.
165
166      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,zmasq,zmasq_2d)
167
168      DO i=1,ip1jmp1-1
169        IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
170                zmasqu(i)=0.
171        ENDIF
172      END DO
173      DO i=iip1,ip1jmp1,iip1
174        zmasqu(i)=zmasqu(i-iim)
175      END DO
176      DO i=1,ip1jm
177        IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN 
178                zmasqv(i)=0.
179        END IF
180      END DO
181
182
183! Coriolis (pour Ekman transp.)
184      eps=1e-5
185!      CALL getin('slab_eps',eps)
186!      print *,'epsilon=',eps
187      ff=fext*unsairez       
188      DO i=1,ip1jm
189        unsev(i)=eps/(ff(i)*ff(i)+eps**2)
190        unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2)
191      ENDDO
192      CALL gr_v_scal(1,unsfv,unsfu)
193      CALL gr_v_scal(1,unsev,unseu)
194! Alpha variable?
195!      alpha_var=.FALSE.
196!      CALL getin('slab_alphav',alpha_var)
197
198
199     
200  END SUBROUTINE init_masquv
201
202
203
204  SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
205 
206      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
207      USE slab_ice_h, ONLY: noceanmx
208
209      IMPLICIT NONE
210     
211      INTEGER,INTENT(IN) :: ngrid
212      INTEGER ij
213      REAL txv((nbp_lon+1)*(nbp_lat-1)),fluxm((nbp_lon+1)*(nbp_lat-1))
214      REAL tyv((nbp_lon+1)*(nbp_lat-1))
215      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),noceanmx)
216      REAL fluxtz((nbp_lon+1)*nbp_lat,noceanmx)
217      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
218      REAL fluxz((nbp_lon+1)*nbp_lat),fluxv((nbp_lon+1)*nbp_lat)
219      REAL dt((nbp_lon+1)*nbp_lat,noceanmx),ts((nbp_lon+1)*nbp_lat,noceanmx)
220      REAL tx_phy(ngrid),ty_phy(ngrid)
221      REAL dt_phy(ngrid,noceanmx),ts_phy(ngrid,noceanmx)
222
223      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
224
225      iim=nbp_lon
226      iip1=nbp_lon+1
227      iip2=nbp_lon+2
228      jjp1=nbp_lat
229      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
230      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
231      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
232
233! Convert taux,y to 2D  scalar grid
234      ! north and south poles tx,ty no meaning
235      tx_phy(1)=0.
236      tx_phy(ngrid)=0.
237      ty_phy(1)=0.
238      ty_phy(ngrid)=0.
239      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,tx_phy,txu)
240      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,ty_phy,tyu)
241         
242! Divide taux,y by f or eps, and convert to 2D u,v grids
243! (Arakawa C grid)
244      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
245      CALL gr_scal_v(1,tyu,tyv)
246      fluxm=tyv*unsev-txv*unsfv
247!      fluxm=-txv*unsfv
248! Zonal flux
249      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
250      CALL gr_scal_u(1,tyu,tyu)
251      fluxz=tyu*unsfu+txu*unseu
252!      fluxz=tyu*unsfu
253           
254! Convert temperature to 2D grid
255      CALL gr_fi_dyn(2,ngrid,iip1,jjp1,ts_phy,ts)
256       
257! Flux de masse
258      fluxm=fluxm*cv*cuvsurcv*zmasqv
259      fluxz=fluxz*cu*cvusurcu*zmasqu
260! Flux de masse vertical
261      DO ij=iip2,ip1jm
262        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
263      ENDDO
264      DO ij=iip1,ip1jmi1,iip1
265         fluxv(ij+1)=fluxv(ij+iip1)
266      END DO
267! Poles
268      fluxv(1)=-SUM(fluxm(1:iim))     
269      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
270      fluxv=fluxv
271
272! Meridional heat fluxes
273      DO ij=1,ip1jm
274          ! centered scheme
275          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
276          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
277      END DO
278
279! Zonal heat fluxes
280! Schema upstream     
281      fluxtz(1:iip1,:)=0 ! no zonal heat flux at north pole
282      DO ij=iip2,ip1jm
283      IF (fluxz(ij).GE.0.) THEN
284             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
285             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
286      ELSE
287             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
288             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
289      ENDIF
290      END DO
291      DO ij=iip1*2,ip1jmp1,iip1
292             fluxtz(ij,:)=fluxtz(ij-iim,:)
293      END DO
294                   
295! Calcul de dT
296      DO ij=iip2,ip1jm
297         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:)   &
298                  +fluxtm(ij,:)-fluxtm(ij-iip1,:)
299         IF (fluxv(ij).GT.0.) THEN
300           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
301           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
302         ELSE
303           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
304           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
305         ENDIF
306         dt(ij,:)=dt(ij,:)/aire(ij)
307      END DO
308      DO ij=iip1,ip1jmi1,iip1
309         dt(ij+1,:)=dt(ij+iip1,:)
310      END DO
311! Pôles
312      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
313        IF (fluxv(1).GT.0.) THEN
314          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
315          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
316        ELSE
317          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
318          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
319        ENDIF
320      dt(1,:)=dt(1,:)/apoln
321      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
322       IF (fluxv(ip1jmp1).GT.0.) THEN
323         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
324         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
325       ELSE
326         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
327         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
328       ENDIF
329      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
330      dt(2:iip1,1)=dt(1,1)
331      dt(2:iip1,2)=dt(1,2)
332      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
333      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
334
335! Retour grille physique
336      CALL gr_dyn_fi(2,iip1,jjp1,ngrid,dt,dt_phy)
337
338
339  END SUBROUTINE slab_ekman2
340
341!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
342
343  SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
344  ! Transfer a variable on global "physics" grid to global "dynamics" grid
345  IMPLICIT NONE
346
347  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
348  REAL,INTENT(IN) :: pfi(ngrid,nfield)
349  REAL,INTENT(OUT) :: pdyn(im,jm,nfield)
350 
351  INTEGER :: i,j,ifield,ig
352 
353  DO ifield=1,nfield
354    ! Handle poles
355    DO i=1,im
356      pdyn(i,1,ifield)=pfi(1,ifield)
357      pdyn(i,jm,ifield)=pfi(ngrid,ifield)
358    ENDDO
359    ! Other points
360    DO j=2,jm-1
361      ig=2+(j-2)*(im-1)
362      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
363      pdyn(im,j,ifield)=pdyn(1,j,ifield)
364    ENDDO
365  ENDDO ! of DO ifield=1,nfield
366 
367  END SUBROUTINE gr_fi_dyn
368 
369
370
371  SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
372  ! Transfer a variable on global "dynamics" grid to global "physics" grid
373  IMPLICIT NONE
374 
375  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
376  REAL,INTENT(IN) :: pdyn(im,jm,nfield)
377  REAL,INTENT(OUT) :: pfi(ngrid,nfield)
378 
379  INTEGER j,ifield,ig
380
381  ! Sanity check:
382  IF(ngrid.NE.2+(jm-2)*(im-1)) THEN
383    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
384    STOP
385  ENDIF
386 
387  ! Handle poles
388  CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
389  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
390  ! Other points
391  DO ifield=1,nfield
392    DO j=2,jm-1
393      ig=2+(j-2)*(im-1)
394      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
395    ENDDO
396  ENDDO
397
398  END SUBROUTINE gr_dyn_fi
399
400
401
402  SUBROUTINE  grad(klevel,pg,pgx,pgy)
403  ! compute the covariant components x,y of the gradient of pg
404  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
405  IMPLICIT NONE
406 
407  INTEGER,INTENT(IN) :: klevel
408  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
409  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
410  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
411
412  INTEGER :: l,ij
413  INTEGER :: iim,iip1,ip1jm,ip1jmp1
414 
415  iim=nbp_lon
416  iip1=nbp_lon+1
417  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
418  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
419 
420  DO l=1,klevel
421    DO ij=1,ip1jmp1-1
422      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
423    ENDDO
424    ! correction for pgx(ip1,j,l) ...
425    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
426    DO ij=iip1,ip1jmp1,iip1
427      pgx(ij,l)=pgx(ij-iim,l)
428    ENDDO
429    DO ij=1,ip1jm
430      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
431    ENDDO
432  ENDDO
433 
434  END SUBROUTINE grad
435
436
437
438  SUBROUTINE diverg(klevel,x,y,div)
439  ! compute the divergence of a vector field of components
440  ! x,y. y and y being covriant components
441  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
442  IMPLICIT NONE
443 
444  INTEGER,INTENT(IN) :: klevel
445  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
446  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
447  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
448 
449  INTEGER :: l,ij
450  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
451 
452  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
453  REAL :: sumypn,sumyps
454  REAL,EXTERNAL :: SSUM
455 
456  iim=nbp_lon
457  iip1=nbp_lon+1
458  iip2=nbp_lon+2
459  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
460  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
461  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
462 
463  DO l=1,klevel
464    DO ij=iip2,ip1jm-1
465      div(ij+1,l)= &
466        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
467        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
468    ENDDO
469    ! correction for div(1,j,l) ...
470    ! ... div(1,j,l)= div(iip1,j,l) ...
471    DO ij=iip2,ip1jm,iip1
472      div(ij,l)=div(ij+iim,l)
473    ENDDO
474    ! at the poles
475    DO ij=1,iim
476      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
477      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
478    ENDDO
479    sumypn=SSUM(iim,aiy1,1)/apoln
480    sumyps=SSUM(iim,aiy2,1)/apols
481    DO ij=1,iip1
482      div(ij,l)=-sumypn
483      div(ij+ip1jm,l)=sumyps
484    ENDDO
485  ENDDO ! of DO l=1,klevel
486 
487  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
488 
489  DO l=1,klevel
490    DO ij=iip2,ip1jm
491      div(ij,l)=div(ij,l)*unsaire(ij)
492    ENDDO
493  ENDDO
494 
495  END SUBROUTINE diverg
496
497
498
499  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
500  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
501  IMPLICIT NONE
502 
503  INTEGER,INTENT(IN) :: nx
504  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
505  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
506 
507  INTEGER :: l,ij
508  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
509 
510  iip1=nbp_lon+1
511  iip2=nbp_lon+2
512  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
513  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
514 
515  DO l=1,nx
516    DO ij=iip2,ip1jm
517      x_scal(ij,l)= &
518                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
519                  /(airev(ij-iip1)+airev(ij))
520    ENDDO
521    DO ij=1,iip1
522      x_scal(ij,l)=0.
523    ENDDO
524    DO ij=ip1jm+1,ip1jmp1
525      x_scal(ij,l)=0.
526    ENDDO
527  ENDDO
528
529  END SUBROUTINE gr_v_scal
530
531  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
532  ! convert values from scalar points to v points on C-grid
533  ! used to compute wind stress at V points
534  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
535  IMPLICIT NONE
536
537  INTEGER,INTENT(IN) :: nx ! number of levels or fields
538  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
539  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
540
541  INTEGER :: l,ij
542  INTEGER :: iip1,ip1jm
543
544  iip1=nbp_lon+1
545  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
546
547      DO l=1,nx
548        DO ij=1,ip1jm
549          x_v(ij,l)= &
550            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
551            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
552            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
553        ENDDO
554      ENDDO
555
556  END SUBROUTINE gr_scal_v
557
558  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
559  ! convert values from scalar points to U points on C-grid
560  ! used to compute wind stress at U points
561  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
562  IMPLICIT NONE
563
564  INTEGER,INTENT(IN) :: nx
565  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
566  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
567
568  INTEGER :: l,ij
569  INTEGER :: iip1,jjp1,ip1jmp1
570
571  iip1=nbp_lon+1
572  jjp1=nbp_lat
573  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
574
575  DO l=1,nx
576     DO ij=1,ip1jmp1-1
577        x_u(ij,l)= &
578         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
579         /(aire(ij)+aire(ij+1))
580     ENDDO
581  ENDDO
582
583  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
584
585  END SUBROUTINE gr_scal_u
586 
587END MODULE surf_heat_transp_mod
588
589
590
Note: See TracBrowser for help on using the repository browser.