source: dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90 @ 4218

Last change on this file since 4218 was 4215, checked in by dubos, 6 years ago

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