Ignore:
Timestamp:
Jan 7, 2020, 5:22:37 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup phyparam (TBC)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90

    r4213 r4215  
    11MODULE phyparam_mod
     2  USE callkeys
     3  USE comgeomfi
    24  IMPLICIT NONE
    35  PRIVATE
    46  SAVE
    57
    6   REAL, PARAMETER :: solarcst=1370., stephan=5.67e-08
     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
    711   
    812  REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
    … …  
    3135       &            pw, &
    3236       &            pdu,pdv,pdt,pdq,pdpsrf)
    33 
    34     USE callkeys
    35     USE comsaison
    36     USE dimphy
    37     USE comgeomfi
    38     USE phys_const
    39     USE planet
    40     USE astronomy
     37    USE phys_const, ONLY : g, rcp, r, unjours
     38    USE surface,    ONLY : soil
    4139    USE turbulence, ONLY : vdif
    4240    USE convection, ONLY : convadj
    43     USE solar, ONLY : solang, zenang, mucorr
    44     USE radiative_sw, ONLY : sw
    45     USE radiative_lw, ONLY : lw
    46     USE surface
    4741
    4842    !=======================================================================
    49     !   Organisation of the physical parametrisations of the LMD
     43    !   Top routine of the physical parametrisations of the LMD
    5044    !   20 parameters GCM for planetary atmospheres.
    5145    !   It includes:
    … …  
    7266         pphi(ngrid,nlayer),    & ! Geopotential at the middle of the layers (m2s-2)
    7367         pphis(ngrid),          & ! surface geopotential (unused)
     68         presnivs(nlayer),      &
    7469         pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
    7570         pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
    … …  
    9590    REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
    9691    REAL zflubid(ngrid),zpmer(ngrid)
    97     REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
     92    REAL zpopsk(ngrid,nlayer)
    9893    REAL zdum1(ngrid,nlayer)
    9994    REAL zdum2(ngrid,nlayer)
    10095    REAL zdum3(ngrid,nlayer)
    101     REAL ztim1,ztim2,ztim3
    102     REAL zls,zinsol
    10396    REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
    10497    REAL zfluxsw(ngrid),zfluxlw(ngrid)
    … …  
    109102    !   ----------------------
    110103   
    111     REAL zps_av
    112    
    113     INTEGER        longcles
    114     PARAMETER    ( longcles = 20 )
    115     REAL clesphy0( longcles      )
    116     REAL presnivs(nlayer)
    117        
    118104    print*,'OK DANS PHYPARAM'
    119105    print*,'nq ',nq
    … …  
    123109         &      size(pdq),size(lati),size(pq),size(factq)
    124110   
    125     nlevel=nlayer+1
    126     igout=ngrid/2+1   
    127     zday=rjourvrai+gmtime
    128 
    129     pdq(:,:,:) = 0. ! we do not use tracers in this physics package
    130 
    131     !-----------------------------------------------------------------------
    132     !    1. Initialisations :
    133     !    --------------------
    134    
    135     IF(firstcall) THEN
    136        
    137        print*,'AKk',ngrid,nsoilmx
    138        allocate(tsurf(ngrid))
    139        print*,'AKa'
    140        allocate (tsoil(ngrid,nsoilmx))
    141        print*,'AKb'
    142        allocate (rnatur(ngrid))
    143        print*,'AK2'
    144        allocate(capcal(ngrid),fluxgrd(ngrid))
    145        print*,'AK3'
    146        allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
    147        print*,'AK4'
    148        allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
    149        print*,'AK5'
    150        allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
    151        print*,'AK6'
    152        
    153        
    154        PRINT*,'FIRSTCALL  '
    155        
    156        !         zday_last=rjourvrai
    157        zday_last=zday-ptimestep/unjours
    158        rnatur=1.
    159        emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
    160        inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
    161        albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
    162        z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
    163        q2=1.e-10
    164        q2l=1.e-10
    165        tsurf=300.
    166        tsoil=300.
    167        print*,tsoil(ngrid/2+1,nsoilmx/2+2)
    168        print*,'TS ',tsurf(igout),tsoil(igout,5)
    169        CALL iniorbit
    170        
    171        if (.not.callrad) then
    172           do ig=1,ngrid
    173              fluxrad(ig)=0.
    174           enddo
    175        endif
    176        
    177        !     print*,'OK PHYPARAM 1 '
    178        IF(callsoil) THEN
    179           CALL soil(ngrid,nsoilmx,firstcall,inertie, &
    180           &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
    181        ELSE
    182           PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
    183           DO ig=1,ngrid
    184              capcal(ig)=1.e5
    185              fluxgrd(ig)=0.
    186           ENDDO
    187        ENDIF
    188        !        CALL inifrict(ptimestep)
    189        icount=0
    190        
    191        PRINT*,'FIRSTCALL B '
    192        print*,'INIIO avant iophys_ini '
    193        call iophys_ini('phys.nc    ',presnivs)
    194        
    195     ENDIF
    196     icount=icount+1
    197    
    198     PRINT*,'FIRSTCALL AP '
    199111    IF (ngrid.NE.ngridmax) THEN
    200112       PRINT*,'STOP in inifis'
    … …  
    204116       STOP
    205117    ENDIF
    206    
    207     DO l=1,nlayer
    208        DO ig=1,ngrid
    209           pdv(ig,l)=0.
    210           pdu(ig,l)=0.
    211           pdt(ig,l)=0.
    212        ENDDO
    213     ENDDO
    214    
    215     DO ig=1,ngrid
    216        pdpsrf(ig)=0.
    217        zflubid(ig)=0.
    218        zdtsrf(ig)=0.
    219     ENDDO
     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.
    220164   
    221165    !-----------------------------------------------------------------------
    … …  
    243187    DO l=1,nlayer
    244188       DO ig=1,ngrid
    245           zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    246        ENDDO
    247     ENDDO
    248     DO l=1,nlayer
    249        DO ig=1,ngrid
     189          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure
    250190          zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    251191       ENDDO
    … …  
    256196    !    ---------------------------------------
    257197   
    258     !        print*,'callrad0'
    259     IF(callrad) THEN
    260        !     print*,'callrad'
    261        
    262        !   WARNING !!! on calcule le ray a chaque appel
    263        !        IF( MOD(icount,iradia).EQ.0) THEN
    264        
    265        CALL solarlong(zday,zls)
    266        CALL orbite(zls,dist_sol,declin)
    267        !      declin=0.
    268        !     print*,'ATTENTIOn : pdeclin = 0','  L_s=',zls
    269        print*,'diurnal=',diurnal
    270        IF(diurnal) THEN
    271           if ( 1.eq.1 ) then
    272              ztim1=SIN(declin)
    273              ztim2=COS(declin)*COS(2.*pi*(zday-.5))
    274              ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
    275              CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
    276              &         ztim1,ztim2,ztim3,                   &
    277              &         mu0,fract)
    278           else
    279              zdtime=ptimestep*float(iradia)
    280              CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
    281              print*,'ZENANG '
    282           endif
    283          
    284           IF(lverbose) THEN
    285              PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
    286              PRINT*,zday, declin,                       &
    287              &            sinlon(igout),coslon(igout),  &
    288              &            sinlat(igout),coslat(igout)
    289           ENDIF
    290        ELSE
    291           print*,'declin,ngrid,rad',declin,ngrid,rad
    292           CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
    293        ENDIF
    294        
    295        !    2.2 Calcul of the radiative tendencies and fluxes:
    296        !    --------------------------------------------------
    297        
    298        !  2.1.2 levels
    299        
    300        zinsol=solarcst/(dist_sol*dist_sol)
    301        zps_av=1.e5
    302        if (firstcall) then
    303           print*,'WARNING ps_rad impose'
    304        endif
    305        CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
    306        &              pplev,zps_av,                 &
    307        &              mu0,fract,zinsol,             &
    308        &              zfluxsw,zdtsw,                &
    309        &              lverbose)
    310        
    311        CALL lw(ngrid,nlayer,coefir,emissiv, &
    312        &             pplev,zps_av,tsurf,pt, &
    313        &             zfluxlw,zdtlw,         &
    314        &             lverbose)
    315 
    316 
    317        !    2.4 total flux and tendencies:
    318        !    ------------------------------
    319        
    320        !    2.4.1 fluxes
    321        
    322        DO ig=1,ngrid
    323           fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
    324                &         +zfluxsw(ig)*(1.-albedo(ig))
    325           zplanck(ig)=tsurf(ig)*tsurf(ig)
    326           zplanck(ig)=emissiv(ig)* &
    327           &         stephan*zplanck(ig)*zplanck(ig)
    328           fluxrad(ig)=fluxrad(ig)-zplanck(ig)
    329        ENDDO
    330        
    331        !    2.4.2 temperature tendencies
    332        
    333        DO l=1,nlayer
    334           DO ig=1,ngrid
    335              dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
    336           ENDDO
    337        ENDDO
    338        
    339        !        ENDIF
    340        
    341        !    2.5 Transformation of the radiative tendencies:
    342        !    -----------------------------------------------
    343        
    344        DO l=1,nlayer
    345           DO ig=1,ngrid
    346              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
    347           ENDDO
    348        ENDDO
    349        
    350        IF(lverbose) THEN
    351           PRINT*,'Diagnotique for the radaition'
    352           PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
    353           PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
    354           &           fract(igout),                       &
    355           &           fluxrad(igout),zplanck(igout)
    356           PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
    357           PRINT*,'unjours',unjours
    358           DO l=1,nlayer
    359              WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
    360              &         pplay(igout,l),pplev(igout,l),  &
    361              &         zdtsw(igout,l),zdtlw(igout,l)
    362           ENDDO
    363        ENDIF
    364        
    365     ENDIF
     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)
    366201
    367202    !-----------------------------------------------------------------------
    … …  
    444279
    445280    !-----------------------------------------------------------------------
    446     !   On incremente les tendances physiques de la temperature du sol:
     281    !   On ajoute les tendances physiques a la temperature du sol:
    447282    !   ---------------------------------------------------------------
    448283   
    … …  
    478313       
    479314       do ig=1,ngridmax
    480           zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*285.))
     315          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
    481316       enddo
    482317       
    … …  
    523358  END SUBROUTINE phyparam
    524359
    525 
    526   SUBROUTINE alloc_phyparam
     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   
    527509  END SUBROUTINE alloc_phyparam
     510
    528511END MODULE phyparam_mod
Note: See TracChangeset for help on using the changeset viewer.