source: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90 @ 4225

Last change on this file since 4225 was 4223, checked in by dubos, 6 years ago

simple_physics : cleanup phyparam and iniphyparam
File size: 15.6 KB
Line 
1MODULE phyparam_mod
2  USE callkeys
3  USE comgeomfi
4  IMPLICIT NONE
5  PRIVATE
6  SAVE
7
8  REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, &
9       ps_rad=1.e5, height_scale=10000., ref_temp=285., &
10       capcal_nosoil=1e5
11   
12  REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
13       capcal(:),fluxgrd(:), &
14       dtrad(:,:),fluxrad(:), &
15       q2(:,:),q2l(:,:), &
16       albedo(:),emissiv(:),z0(:),inertie(:)
17  !$OMP THREADPRIVATE( tsurf,tsoil,rnatur)
18  !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
19  !$OMP THREADPRIVATE( q2,q2l)
20  !$OMP THREADPRIVATE( albedo,emissiv,z0,inertie)         
21
22  INTEGER :: icount
23  REAL    :: zday_last
24  !$OMP THREADPRIVATE( icount,zday_last)
25
26  PUBLIC :: phyparam
27
28CONTAINS
29
30  SUBROUTINE phyparam(ngrid,nlayer,             &
31       &            firstcall,lastcall,         &
32       &            rjourvrai,gmtime,ptimestep, &
33       &            pplev,pplay,pphi,           &
34       &            pu,pv,pt,                   &
35       &            pdu,pdv,pdt,pdpsrf)
36    USE phys_const, ONLY : g, rcp, r, unjours
37    USE surface,    ONLY : soil
38    USE turbulence, ONLY : vdif
39    USE convection, ONLY : convadj
40    USE writefield_mod, ONLY : writefield
41
42    !=======================================================================
43    !   Top routine of the physical parametrisations of the LMD
44    !   20 parameters GCM for planetary atmospheres.
45    !   It includes:
46    !   raditive transfer (long and shortwave) for CO2 and dust.
47    !   vertical turbulent mixing
48    !   convective adjsutment
49    !
50    !   author: Frederic Hourdin 15 / 10 /93
51    !=======================================================================
52
53    INTEGER, INTENT(IN) ::      &
54         ngrid,                 & ! Size of the horizontal grid.
55         nlayer                   ! Number of vertical layers.
56    LOGICAL, INTENT(IN) ::      &
57         firstcall,             & ! True at the first call
58         lastcall                 ! True at the last call
59    REAL, INTENT(IN)    ::      &
60         rjourvrai,             & ! Number of days counted from the North. Spring equinox
61         gmtime,                & ! time of the day in seconds
62         ptimestep,             & ! timestep (s)
63         pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa)
64         pplay(ngrid,nlayer),   & ! Pressure at the middle of the layers (Pa)
65         pphi(ngrid,nlayer),    & ! Geopotential at the middle of the layers (m2s-2)
66         pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
67         pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
68         pt(ngrid,nlayer)         ! Temperature (K)
69    REAL, INTENT(OUT)   ::      & ! output : physical tendencies
70         pdu(ngrid,nlayer),     &
71         pdv(ngrid,nlayer),     &
72         pdt(ngrid,nlayer),     &
73         pdpsrf(ngrid)
74   
75    !    Local variables :
76    REAL, DIMENSION(ngrid) :: mu0,fract
77    INTEGER :: j,l,ig,nlevel,igout
78    !
79    REAL :: zday, zdtime
80    REAL zh(ngrid,nlayer),z1,z2
81    REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
82    REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
83    REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
84    REAL zflubid(ngrid),zpmer(ngrid)
85    REAL zpopsk(ngrid,nlayer)
86    REAL zdum1(ngrid,nlayer)
87    REAL zdum2(ngrid,nlayer)
88    REAL zdum3(ngrid,nlayer)
89    REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
90    REAL zfluxsw(ngrid),zfluxlw(ngrid)
91   
92    print*,'OK DANS PHYPARAM'
93    print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
94    print*,'nlayer',nlayer
95   
96    IF (ngrid.NE.ngridmax) THEN
97       PRINT*,'STOP in inifis'
98       PRINT*,'Probleme de dimenesions :'
99       PRINT*,'ngrid     = ',ngrid
100       PRINT*,'ngridmax   = ',ngridmax
101       STOP
102    ENDIF
103
104    nlevel=nlayer+1
105    igout=ngrid/2+1   
106    zday=rjourvrai+gmtime
107
108    !-----------------------------------------------------------------------
109    !    0. Allocate and initialize at first call
110    !    --------------------
111   
112    IF(firstcall) THEN
113       !         zday_last=rjourvrai
114       zday_last=zday-ptimestep/unjours
115       CALL alloc_phyparam(ngrid, nlayer, igout)
116
117       !     print*,'OK PHYPARAM 1 '
118       IF(callsoil) THEN
119          CALL soil(ngrid,nsoilmx,firstcall,inertie, &
120               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
121          ! NB : this call to soil also performs some calculations, see surface.F90
122       ELSE
123          PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
124          DO ig=1,ngrid
125             capcal(ig)  = capcal_nosoil
126             fluxgrd(ig) = 0.
127          ENDDO
128       ENDIF
129    ENDIF
130
131    !-----------------------------------------------------------------------
132    !    1. Initialisations :
133    !    --------------------
134   
135    icount=icount+1
136   
137    pdv(:,:)  = 0.
138    pdu(:,:)  = 0.
139    pdt(:,:)  = 0.
140    pdpsrf(:) = 0.
141    zflubid(:)= 0.
142    zdtsrf(:) = 0.
143   
144    !-----------------------------------------------------------------------
145    !   calcul du geopotentiel aux niveaux intercouches
146    !   ponderation des altitudes au niveau des couches en dp/p
147   
148    DO l=1,nlayer
149       DO ig=1,ngrid
150          zzlay(ig,l)=pphi(ig,l)/g
151       ENDDO
152    ENDDO
153    DO ig=1,ngrid
154       zzlev(ig,1)=0.
155    ENDDO
156    DO l=2,nlayer
157       DO ig=1,ngrid
158          z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
159          z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
160          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
161       ENDDO
162    ENDDO
163   
164    !-----------------------------------------------------------------------
165    !   Transformation de la temperature en temperature potentielle
166    DO l=1,nlayer
167       DO ig=1,ngrid
168          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure
169          zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
170       ENDDO
171    ENDDO
172   
173    !-----------------------------------------------------------------------
174    !    2. Calcul of the radiative tendencies :
175    !    ---------------------------------------
176   
177    IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
178         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
179         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
180
181    !-----------------------------------------------------------------------
182    !    3. Vertical diffusion (turbulent mixing):
183    !    -----------------------------------------
184    !
185    IF(calldifv) THEN
186       
187       DO ig=1,ngrid
188          zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
189       ENDDO
190       
191       zdum1(:,:)=0.
192       zdum2(:,:)=0.
193       
194       do l=1,nlayer
195          do ig=1,ngrid
196             zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
197          enddo
198       enddo
199       
200       CALL vdif(ngrid,nlayer,zday,                 &
201       &        ptimestep,capcal,z0,                &
202       &        pplay,pplev,zzlay,zzlev,            &
203       &        pu,pv,zh,tsurf,emissiv,             &
204       &        zdum1,zdum2,zdum3,zflubid,          &
205       &        zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,   &
206       &        lverbose)
207
208       DO l=1,nlayer
209          DO ig=1,ngrid
210             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
211             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
212             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
213          ENDDO
214       ENDDO
215       
216       DO ig=1,ngrid
217          zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
218       ENDDO
219       
220    ELSE
221       DO ig=1,ngrid
222          zdtsrf(ig)=zdtsrf(ig)+ &
223          &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
224       ENDDO
225    ENDIF
226    !
227    !-----------------------------------------------------------------------
228    !   4. Dry convective adjustment:
229    !   -----------------------------
230   
231    IF(calladj) THEN
232       
233       DO l=1,nlayer
234          DO ig=1,ngrid
235             zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
236          ENDDO
237       ENDDO
238       
239       zdufr(:,:)=0.
240       zdvfr(:,:)=0.
241       zdhfr(:,:)=0.
242       
243       CALL convadj(ngrid,nlayer,ptimestep, &
244       &                pplay,pplev,zpopsk, &
245       &                pu,pv,zh,           &
246       &                pdu,pdv,zdum1,      &
247       &                zdufr,zdvfr,zdhfr)
248       
249       DO l=1,nlayer
250          DO ig=1,ngrid
251             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
252             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
253             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
254          ENDDO
255       ENDDO
256       
257    ENDIF
258
259    !-----------------------------------------------------------------------
260    !   On ajoute les tendances physiques a la temperature du sol:
261    !   ---------------------------------------------------------------
262   
263    DO ig=1,ngrid
264       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
265    ENDDO
266   
267    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
268   
269    !-----------------------------------------------------------------------
270    !   soil temperatures:
271    !   --------------------
272   
273    IF (callsoil) THEN
274       CALL soil(ngrid,nsoilmx,.false.,inertie, &
275       &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
276       IF(lverbose) THEN
277          PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
278          PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), &
279               &          zdtsrf(igout),ptimestep
280       ENDIF
281    ENDIF
282   
283    !-----------------------------------------------------------------------
284    !   sorties:
285    !   --------
286   
287    print*,'zday, zday_last ',zday,zday_last,icount
288    if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
289       print*,'zday, zday_last SORTIE ',zday,zday_last
290       zday_last=zday
291       !  Ecriture/extension de la coordonnee temps
292       
293       do ig=1,ngridmax
294          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
295       enddo
296       
297       call writefield('u','Vent zonal moy','m/s',pu)
298       call writefield('v','Vent meridien moy','m/s',pv)
299       call writefield('temp','Temperature','K',pt)
300       call writefield('geop','Geopotential','m2/s2',pphi)
301       call writefield('plev','plev','Pa',pplev(:,1:nlayer))
302       
303       call writefield('du','du',' ',pdu)
304       call writefield('dv','du',' ',pdv)
305       call writefield('dt','du',' ',pdt)
306       call writefield('dtsw','dtsw',' ',zdtsw)
307       call writefield('dtlw','dtlw',' ',zdtlw)
308       
309       call writefield('ts','Surface temper','K',tsurf)
310       call writefield('coslon','coslon',' ',coslon)
311       call writefield('sinlon','sinlon',' ',sinlon)
312       call writefield('coslat','coslat',' ',coslat)
313       call writefield('sinlat','sinlat',' ',sinlat)
314       call writefield('mu0','mu0',' ',mu0)
315       call writefield('alb','alb',' ',albedo)
316       call writefield('fract','fract',' ',fract)
317       call writefield('ps','Surface pressure','Pa',pplev(:,1))
318       call writefield('slp','Sea level pressure','Pa',zpmer)
319       call writefield('swsurf','SW surf','Pa',zfluxsw)
320       call writefield('lwsurf','LW surf','Pa',zfluxlw)
321       
322    endif
323   
324  END SUBROUTINE phyparam
325
326  SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, &
327       gmtime, zdtime, zday, pplev, pplay, pt, &
328       pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
329    USE planet
330    USE phys_const,   ONLY : planet_rad, unjours
331    USE astronomy,    ONLY : orbite, solarlong
332    USE solar,        ONLY : solang, zenang, mucorr
333    USE radiative_sw, ONLY : sw
334    USE radiative_lw, ONLY : lw
335
336    INTEGER, INTENT(IN) :: ngrid, igout, nlayer
337    REAL, INTENT(IN)    :: gmtime, zdtime, zday, &
338         &                 pplev(ngrid,nlayer+1), pplay(ngrid, nlayer), &
339         &                 pt(ngrid, nlayer+1)
340    REAL, INTENT(INOUT) :: pdt(ngrid,nlayer)
341    REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), zfluxlw(ngrid), &
342         zdtsw(ngrid,nlayer), zfluxsw(ngrid), mu0(ngrid)
343    REAL, DIMENSION(ngrid) :: fract
344    REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3, dist_sol, declin
345    REAL :: zplanck(ngrid)
346    INTEGER :: ig, l
347
348    !    2.1 Insolation
349    !    --------------------------------------------------
350   
351    CALL solarlong(zday,zls)
352    CALL orbite(zls,dist_sol,declin)
353   
354    IF(diurnal) THEN
355       IF ( .TRUE. ) then
356          ztim1=SIN(declin)
357          ztim2=COS(declin)*COS(2.*pi*(zday-.5))
358          ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
359          CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
360               &         ztim1,ztim2,ztim3,                   &
361               &         mu0,fract)
362       ELSE
363          CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
364          print*,'ZENANG '
365       ENDIF
366       
367       IF(lverbose) THEN
368          PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
369          PRINT*,zday, declin,                       &
370               &            sinlon(igout),coslon(igout),  &
371               &            sinlat(igout),coslat(igout)
372       ENDIF
373    ELSE
374       print*,'declin,ngrid,planet_rad',declin,ngrid,planet_rad
375       CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
376    ENDIF
377   
378    zinsol=solarcst/(dist_sol*dist_sol)
379   
380    !    2.2 Radiative tendencies and fluxes:
381    !    --------------------------------------------------
382   
383    CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
384         &              pplev,ps_rad,                 &
385         &              mu0,fract,zinsol,             &
386         &              zfluxsw,zdtsw,                &
387         &              lverbose)
388   
389    CALL lw(ngrid,nlayer,coefir,emissiv, &
390         &             pplev,ps_rad,tsurf,pt, &
391         &             zfluxlw,zdtlw,         &
392         &             lverbose)
393   
394    !    2.4 surface fluxes
395    !    ------------------------------
396   
397    DO ig=1,ngrid
398       fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
399            &         +zfluxsw(ig)*(1.-albedo(ig))
400       tsurf2 = tsurf(ig)*tsurf(ig)
401       zplanck(ig)=emissiv(ig)*stephan*tsurf2*tsurf2
402       fluxrad(ig)=fluxrad(ig)-zplanck(ig)
403    ENDDO
404   
405    !    2.5 Temperature tendencies
406    !    --------------------------
407   
408    DO l=1,nlayer
409       DO ig=1,ngrid
410          dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
411          pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
412       ENDDO
413    ENDDO
414   
415    IF(lverbose) THEN
416       PRINT*,'Diagnostics for radiation'
417       PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
418       PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
419            &           fract(igout),                       &
420            &           fluxrad(igout),zplanck(igout)
421       PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
422       PRINT*,'unjours',unjours
423       DO l=1,nlayer
424          WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
425               &         pplay(igout,l),pplev(igout,l),  &
426               &         zdtsw(igout,l),zdtlw(igout,l)
427       ENDDO
428    ENDIF
429   
430  END SUBROUTINE radiative_tendencies
431
432  SUBROUTINE alloc_phyparam(ngrid, nlayer, igout)
433    USE surface
434    USE astronomy, ONLY : iniorbit
435    INTEGER, INTENT(IN) :: ngrid, nlayer, igout
436    LOGICAL, PARAMETER :: firstcall=.TRUE.
437
438    print*,'AKk',ngrid,nsoilmx
439    allocate(tsurf(ngrid))
440    print*,'AKa'
441    allocate (tsoil(ngrid,nsoilmx))
442    print*,'AKb'
443    allocate (rnatur(ngrid))
444    print*,'AK2'
445    allocate(capcal(ngrid),fluxgrd(ngrid))
446    print*,'AK3'
447    allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
448    print*,'AK4'
449    allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
450    print*,'AK5'
451    allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
452    print*,'AK6'
453   
454   
455    PRINT*,'FIRSTCALL  '
456   
457    rnatur=1.
458    emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
459    inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
460    albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
461    z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
462    q2=1.e-10
463    q2l=1.e-10
464    tsurf=300.
465    tsoil=300.
466    print*,tsoil(igout,nsoilmx/2+2)
467    print*,'TS ',tsurf(igout),tsoil(igout,5)
468    CALL iniorbit
469   
470    if (.not.callrad) fluxrad(:)=0.
471   
472    icount=0
473   
474  END SUBROUTINE alloc_phyparam
475
476END MODULE phyparam_mod
Note: See TracBrowser for help on using the repository browser.