Changeset 3672 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Mar 6, 2025, 11:25:14 AM (3 months ago)
Author:
afalco
Message:

generic/pluto: tsurf/tsoil taken from phys_state_var_mod rather than declared locally so they are updated in the outputs (profile.out, restart1D.nc).
AF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F

    r3663 r3672  
    3636      use mod_interface_dyn_phys, only: init_interface_dyn_phys
    3737      use inifis_mod, only: inifis
    38       use phys_state_var_mod, only: phys_state_var_init
     38      use phys_state_var_mod, only: phys_state_var_init, tsurf, tsoil
    3939      use physiq_mod, only: physiq
    4040      use restart1D_mod, only: writerestart1D
    … …  
    4848
    4949!==================================================================
    50 !     
     50!
    5151!     Purpose
    5252!     -------
    5353!     Run the physics package of the universal model in a 1D column.
    54 !     
     54!
    5555!     It can be compiled with a command like (e.g. for 25 layers):
    5656!                                  "makegcm -p std -d 25 rcm1d"
    … …  
    9595      REAL play(llm)        ! Pressure at the middle of the layers (Pa)
    9696      REAL plev(llm+1)      ! intermediate pressure levels (pa)
    97       REAL psurf,tsurf(1)     
     97      REAL psurf
    9898      REAL u(llm),v(llm)    ! zonal, meridional wind
    9999      REAL gru,grv          ! prescribed "geostrophic" background wind
    … …  
    101101      REAL,ALLOCATABLE :: q(:,:)      ! tracer mixing ratio (e.g. kg/kg)
    102102      REAL,ALLOCATABLE :: qsurf(:)    ! tracer surface budget (e.g. kg.m-2)
    103       REAL,ALLOCATABLE :: tsoil(:)    ! subsurface soil temperature (K)
    104103!      REAL co2ice               ! co2ice layer (kg.m-2) !not used anymore
    105104      integer :: i_co2_ice=0     ! tracer index of co2 ice
    … …  
    114113      REAL du(llm),dv(llm),dtemp(llm)
    115114      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
    116       REAL dpsurf(1)   
     115      REAL dpsurf(1)
    117116      REAL,ALLOCATABLE :: dq(:,:)
    118117      REAL,ALLOCATABLE :: dqdyn(:,:)
    … …  
    126125      REAL tmp1(0:llm),tmp2(0:llm)
    127126      integer :: nq !=1 ! number of tracers
    128  
     127
    129128      character*2 str2
    130129      character (len=7) :: str7
    … …  
    181180      endif
    182181
    183 ! check if 'rcm1d.def' file is around 
     182! check if 'rcm1d.def' file is around
    184183      open(90,file='rcm1d.def',status='old',form='formatted',
    185184     &     iostat=ierr)
    … …  
    257256      endif
    258257      close(90)
    259      
     258
    260259      ! Initialize dimphy module
    261       call init_dimphy(1,llm) 
    262      
     260      call init_dimphy(1,llm)
     261
    263262      ! now initialize arrays using phys_state_var_init
    264263      ! but first initialise naerkind (from callphys.def)
    265264      naerkind=0 !default
    266265      call getin("naerkind",naerkind)
    267      
     266
    268267      call phys_state_var_init(nq)
    269      
     268
    270269      saveprofile=.false.
    271270      saveprofile=.true.
    … …  
    277276      pi=2.E+0*asin(1.E+0)
    278277
    279 c     Parametres Couche limite et Turbulence 
     278c     Parametres Couche limite et Turbulence
    280279c     --------------------------------------
    281       z0 =  1.e-2                ! surface roughness (m) ~0.01 
    282  
     280      z0 =  1.e-2                ! surface roughness (m) ~0.01
     281
    283282c     propriete optiques des calottes et emissivite du sol
    284283c     ----------------------------------------------------
    285284      emissiv= 0.95              ! Emissivite du sol martien ~.95
    286285      emisice(1)=0.95            ! Emissivite calotte nord
    287       emisice(2)=0.95            ! Emissivite calotte sud 
     286      emisice(2)=0.95            ! Emissivite calotte sud
    288287
    289288      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
    … …  
    294293
    295294c ------------------------------------------------------
    296 c  Load parameters from "run.def" and "gases.def" 
     295c  Load parameters from "run.def" and "gases.def"
    297296c ------------------------------------------------------
    298297
    … …  
    355354            stop
    356355          endif
    357        
     356
    358357          do iq=1,nq
    359358          ! minimal version, just read in the tracer names, 1 per line
    … …  
    391390        nq=0
    392391        ! still, make allocations for 1 dummy tracer
    393         allocate(tname(1)) 
     392        allocate(tname(1))
    394393        allocate(qsurf(1))
    395394        allocate(q(llm,1))
    396395        allocate(dq(llm,1))
    397      
     396
    398397       ! Check that tracer boolean is compliant with number of tracers
    399398       ! -- otherwise there is an error (and more generally we have to be consistent)
    … …  
    413412     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
    414413      phisfi(1)=0.E+0
    415      !!! LATITUDE. can be set. 
     414     !!! LATITUDE. can be set.
    416415      latitude=0 ! default value for latitude
    417416      PRINT *,'latitude (in degrees) ?'
    … …  
    468467          PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
    469468          STOP
    470       ELSE 
     469      ELSE
    471470          PRINT *,"--> year_day = ",year_day
    472471      ENDIF
    … …  
    501500          PRINT *,"STOP. peri_day.gt.year_day"
    502501          STOP
    503       ELSE 
     502      ELSE
    504503          PRINT *,"--> peri_day = ", peri_day
    505       ENDIF 
    506 
    507       obliquit = -99999. 
     504      ENDIF
     505
     506      obliquit = -99999.
    508507      PRINT *,'OBLIQUITY in deg ?'
    509508      call getin("obliquit",obliquit)
    … …  
    513512      ELSE
    514513          PRINT *,"--> obliquit = ",obliquit
    515       ENDIF 
     514      ENDIF
    516515
    517516      if (restart) then
    … …  
    552551c    Date (en sols depuis le solstice de printemps) du debut du run
    553552      !if (restart) then
    554       !  ! day   
     553      !  ! day
    555554      !  ierr=NF90_INQ_VARID(nid_restart1D,'day',did)
    556555      !  if (ierr==NF90_NOERR) then
    … …  
    609608      nday=ndt
    610609
    611       ndt=ndt*day_step     
    612       dtphys=daysec/day_step 
     610      ndt=ndt*day_step
     611      dtphys=daysec/day_step
    613612      write(*,*)"-------------------------------------"
    614613      write(*,*)"-------------------------------------"
    … …  
    635634!!! - read callphys.def
    636635!!! - calculate sine and cosine of longitude and latitude
    637 !!! - calculate mugaz and cp 
     636!!! - calculate mugaz and cp
    638637!!! NB: some operations are being done dummily in inifis in 1D
    639638!!! - physical constants: nevermind, things are done allright below
    … …  
    645644
    646645      nsoil=nsoilmx
    647       allocate(tsoil(nsoilmx))
    648       !! those are defined in comsoil_h.F90
    649       IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
    650       IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
    651       IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
    652646
    653647! At this point, both getin() and getin_p() functions have been used,
    … …  
    698692        ENDDO
    699693      ENDDO
    700      
     694
    701695      DO iq=1,nq
    702696        qsurf(iq) = 0.
    703697      ENDDO
    704      
    705       if (tracer) then ! Initialize tracers here. 
    706              
     698
     699      if (tracer) then ! Initialize tracers here.
     700
    707701         write(*,*) "rcm1d : initializing tracers profiles"
    708702
    … …  
    738732              endif
    739733            else ! restart
    740          
     734
    741735              txt=""
    742736              write(txt,"(a)") tname(iq)
    743737              write(*,*)"  tracer:",trim(txt)
    744              
     738
    745739              ! CO2
    746740              if (txt.eq."co2_ice") then
    747741                 q(:,iq)=0.   ! kg/kg of atmosphere
    748                  qsurf(iq)=0. ! kg/m2 at the surface               
     742                 qsurf(iq)=0. ! kg/m2 at the surface
    749743                 ! Look for a "profile_co2_ice" input file
    750744                 open(91,file='profile_co2_ice',status='old',
    … …  
    760754                 close(91)
    761755              endif ! of if (txt.eq."co2")
    762          
     756
    763757              ! WATER VAPOUR
    764758              if (txt.eq."h2o_vap") then
    765759                 q(:,iq)=0.   ! kg/kg of atmosphere
    766760                 qsurf(iq)=0. ! kg/m2 at the surface
    767                  ! Look for a "profile_h2o_vap" input file   
     761                 ! Look for a "profile_h2o_vap" input file
    768762                 open(91,file='profile_h2o_vap',status='old',
    769763     &           form='formatted',iostat=ierr)
    … …  
    778772                 close(91)
    779773              endif ! of if (txt.eq."h2o_vap")
    780            
     774
    781775              ! WATER ICE
    782776              if (txt.eq."h2o_ice") then
    … …  
    807801                    open(91,file=tracer_profile_file_name,status='old',
    808802     &              form="formatted",iostat=ierr)
    809                     if (ierr .eq. 0) then 
     803                    if (ierr .eq. 0) then
    810804                          read(91,*)qsurf(iq)
    811                           do ilayer=1,nlayer 
     805                          do ilayer=1,nlayer
    812806                                read(91,*)q(ilayer,iq)
    813                           enddo 
    814                     else 
     807                          enddo
     808                    else
    815809                          write(*,*) "No initial profile "
    816810                          write(*,*) " for this tracer :"
    … …  
    818812                    endif
    819813                    close(91)
    820               endif ! (txt .eq. "_vap") 
    821               !_ice 
    822               if((txt.ne."h2o_ice") .and. 
     814              endif ! (txt .eq. "_vap")
     815              !_ice
     816              if((txt.ne."h2o_ice") .and.
    823817     &                        (index(txt,'_ice'   ) /= 0)) then
    824818                    q(:,iq)=0. !kg/kg of atmosphere
    … …  
    827821            endif ! restart
    828822           enddo ! of do iq=1,nq
    829        
     823
    830824      endif ! of tracer
    831825
    … …  
    833827c   ------------------------------------------------------
    834828      ptif=2.E+0*omeg*sinlat(1)
    835  
     829
    836830
    837831c    vent geostrophique
    … …  
    941935
    942936      if(autozlevs)then
    943            
     937
    944938         open(91,file="z2sig.def",form='formatted')
    945939         read(91,*) Hscale
    … …  
    948942         enddo
    949943         close(91)
    950  
    951            
     944
     945
    952946         print*,'Hmax = ',Hmax,' km'
    953947         print*,'Auto-shifting Hscale to:'
    … …  
    955949         Hscale = Hmax / log(psurf/pceil)
    956950         print*,'Hscale = ',Hscale,' km'
    957          
     951
    958952! none of this matters if we dont care about zlay
    959          
     953
    960954      endif
    961955
    … …  
    981975            plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
    982976         ENDDO
    983          
     977
    984978         DO ilayer=1,nlayer
    985979            play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
    986980         ENDDO
    987          
     981
    988982
    989983
    … …  
    10611055      DO isoil=1,nsoil
    10621056         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
    1063          tsoil(isoil)=tsurf(1)  ! soil temperature
     1057         tsoil(1,isoil)=tsurf(1)  ! soil temperature
    10641058      ENDDO
    10651059
    … …  
    11391133     &                  dtphys,time,
    11401134     &                  tsurf,tsoil,emis,albedo,q2,qsurf,
    1141      &                  cloudfrac,totcloudfrac,hice, 
     1135     &                  cloudfrac,totcloudfrac,hice,
    11421136     &                  rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice)
    11431137      endif
    11441138c=======================================================================
    1145 c  BOUCLE TEMPORELLE DU MODELE 1D 
     1139c  BOUCLE TEMPORELLE DU MODELE 1D
    11461140c=======================================================================
    11471141
    … …  
    11621156        ENDIF
    11631157
    1164 c    calcul du geopotentiel 
     1158c    calcul du geopotentiel
    11651159c     ~~~~~~~~~~~~~~~~~~~~~
    11661160
    … …  
    11971191     ,     day,time,dtphys,
    11981192     ,     plev,play,phi,
    1199      ,     u, v,temp, q, 
     1193     ,     u, v,temp, q,
    12001194     ,     w,
    12011195C - sorties
    … …  
    12051199c       evolution du vent : modele 1D
    12061200c       -----------------------------
    1207  
     1201
    12081202c       la physique calcule les derivees temporelles de u et v.
    12091203c       on y rajoute betement un effet Coriolis.
    … …  
    12211215          ENDDO
    12221216c       end if
    1223 c     
     1217c
    12241218c
    12251219c       Calcul du temps au pas de temps suivant
    … …  
    12811275      if (lastcall) then
    12821276        call writerestart1D('restart1D.nc',nlayer,nsoil,day,time,psurf,
    1283      &                                     temp,tsoil,u,v,nq,q)   
     1277     &                                     temp,tsoil,u,v,nq,q)
    12841278      endif
    12851279
    … …  
    12901284c    ========================================================
    12911285      end                       !rcm1d
    1292  
    1293 
     1286
     1287
Note: See TracChangeset for help on using the changeset viewer.