Changeset 2281


Ignore:
Timestamp:
Apr 8, 2020, 6:35:50 PM (5 years ago)
Author:
adelavois
Message:

Mars GCM:
Martian physics is now able to start without startfi.nc
Major update for phyetat0_mod and physiq_mod based on what have been done for the Generic physics

Location:
trunk/LMDZ.MARS
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2280 r2281  
    29322932This program can be used to interpolate GCM files at one terminator (morning or evening) all around the globe.
    29332933+ a fix of localtime.F90 which didn't work for stats files.
     2934
     2935== 08/04/2020 == AD
     2936Martian physics is now able to start without startfi.nc
     2937Major update for phyetat0_mod and physiq_mod based on what have been done for the Generic physics

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r2221 r2281  
    1616     &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
    1717     &   ,satindexco2,rdstorm,slpwind,calllott_nonoro                   &
    18      &   ,latentheat_surfwater,gwd_convective_source
    19      
     18     &   ,latentheat_surfwater,gwd_convective_source,startphy_file
     19
    2020      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    2121     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear,dustinjection
    … …  
    3333
    3434      COMMON/aeroutput/dustiropacity
     35
     36      logical startphy_file
    3537
    3638      logical callemis
    … …  
    8183      integer nircorr
    8284
     85
    8386      character(len=100) dustiropacity
    8487      real               dustrefir

  • trunk/LMDZ.MARS/libf/phymars/conf_phys.F

    r2218 r2281  
    3535! to use  'getin'
    3636      USE ioipsl_getincom, only : getin
     37      USE ioipsl_getin_p_mod, ONLY : getin_p
    3738      use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
    3839     &                       nuice_ref,nuiceco2_ref
    … …  
    9091         write(*,*) " datadir = ",trim(datadir)
    9192
     93         write(*,*) "Initialize physics with startfi.nc file ?"
     94         startphy_file=.true.
     95         call getin_p("startphy_file",startphy_file)
     96         write(*,*) "startphy_file", startphy_file
     97         
    9298         write(*,*) "Run with or without tracer transport ?"
    9399         tracer=.false. ! default value

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r2265 r2281  
    1616                     get_field, get_var, inquire_field, &
    1717                     inquire_dimension, inquire_dimension_length
    18   use ioipsl_getincom, only : getin
    1918  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
    2019  use compute_dtau_mod, only: dtau
     20
     21  USE ioipsl_getin_p_mod, ONLY : getin_p
    2122
    2223  implicit none
    … …  
    3233!         June 2013 TN : Possibility to read files with a time axis
    3334!         November 2013 EM : Enabeling parallel, using iostart module
     35!         March 2020 AD: Enabling initialization of physics without startfi
     36!                        flag: startphy_file
    3437!======================================================================
    3538  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
    … …  
    3942!  ---------
    4043!  inputs:
     44!  logical,intent(in) :: startphy_file ! .true. if reading start file
    4145  character*(*),intent(in) :: fichnom ! "startfi.nc" file
    4246  integer,intent(in) :: tab0
    … …  
    97101
    98102      REAL :: timestart ! to pick which initial state to start from
    99 
    100 ! open physics initial state file:
    101 call open_startphy(fichnom)
    102 
    103 
    104 ! possibility to modify tab_cntrl in tabfi
    105 write(*,*)
    106 write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
    107 call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
    108             p_omeg,p_g,p_mugaz,p_daysec,time0)
    109 
    110 
    111 ! Load surface geopotential:
    112 call get_field("phisfi",phisfi,found)
    113 if (.not.found) then
    114   write(*,*) "phyetat0: Failed loading <phisfi>"
    115   call abort
    116 else
    117   write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
    118              minval(phisfi), maxval(phisfi)
     103      REAL :: surfemis  ! constant emissivity when no startfi
     104      REAL :: surfalbedo  ! constant emissivity when no startfi
     105      CHARACTER(len=5) :: modname="phyetat0"
     106
     107write(*,*) "phyetat0: startphy_file", startphy_file
     108
     109if(.not.startphy_file) then
     110   surfemis = -999
     111   surfalbedo = -999
     112   call getin_p("surfemis",surfemis)
     113   call getin_p("surfalbedo",surfalbedo)
     114   if (surfemis.eq.-999) then
     115      call abort_physic(modname, &
     116               "Missing value for surfemis in def files!",1)
     117   endif
     118   if (surfalbedo.eq.-999) then
     119      call abort_physic(modname, &
     120               "Missing value for surfalbedo in def files!",1)
     121   endif
     122   write(*,*) "phyetat0: surfemis: ", surfemis
     123   write(*,*) "phyetat0: surfalbedo: ", surfalbedo
    119124endif
    120125
    121 
    122 ! Load bare ground albedo:
    123 call get_field("albedodat",albedodat,found)
    124 if (.not.found) then
    125   write(*,*) "phyetat0: Failed loading <albedodat>"
    126   call abort
    127 else
    128   write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
     126if (startphy_file) then
     127   ! open physics initial state file:
     128   call open_startphy(fichnom)
     129   ! possibility to modify tab_cntrl in tabfi
     130   write(*,*)
     131   write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
     132   call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
     133               p_omeg,p_g,p_mugaz,p_daysec,time0)
     134else ! "academic" initialization of planetary parameters via tabfi
     135   call tabfi (0,0,0,day_ini,lmax,p_rad, &
     136               p_omeg,p_g,p_mugaz,p_daysec,time0)
     137endif ! of if (startphy_file)
     138
     139if (startphy_file) then
     140   ! Load surface geopotential:
     141   call get_field("phisfi",phisfi,found)
     142   if (.not.found) then
     143     call abort_physic(modname, &
     144                "phyetat0: Failed loading <phisfi>",1)
     145   endif
     146else
     147  phisfi(:)=0.
     148endif ! of if (startphy_file)
     149write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
     150               minval(phisfi), maxval(phisfi)
     151
     152
     153if (startphy_file) then
     154   ! Load bare ground albedo:
     155   call get_field("albedodat",albedodat,found)
     156   if (.not.found) then
     157     call abort_physic(modname, &
     158                "phyetat0: Failed loading <albedodat>",1)
     159   endif
     160else ! If no startfi file, use parameter surfalbedo in def file
     161  surfalbedo=0.5
     162  call getin_p("surfalbedo",surfalbedo)
     163  print*,"surfalbedo",surfalbedo
     164  albedodat(:)=surfalbedo
     165endif ! of if (startphy_file)
     166write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
    129167             minval(albedodat), maxval(albedodat)
    130 endif
    131168
    132169! ZMEA
    133 call get_field("ZMEA",zmea,found)
    134 if (.not.found) then
    135   write(*,*) "phyetat0: Failed loading <ZMEA>"
    136   call abort
    137 else
    138   write(*,*) "phyetat0: <ZMEA> range:", &
    139              minval(zmea), maxval(zmea)
    140 endif
     170if (startphy_file) then
     171   call get_field("ZMEA",zmea,found)
     172   if (.not.found) then
     173     call abort_physic(modname, &
     174                "phyetat0: Failed loading <ZMEA>",1)
     175   endif
     176else
     177  zmea(:)=0.
     178endif ! of if (startphy_file)
     179write(*,*) "phyetat0: <ZMEA> range:", &
     180            minval(zmea), maxval(zmea)
    141181
    142182
    143183! ZSTD
    144 call get_field("ZSTD",zstd,found)
    145 if (.not.found) then
    146   write(*,*) "phyetat0: Failed loading <ZSTD>"
    147   call abort
    148 else
    149   write(*,*) "phyetat0: <ZSTD> range:", &
    150              minval(zstd), maxval(zstd)
    151 endif
     184if (startphy_file) then
     185   call get_field("ZSTD",zstd,found)
     186   if (.not.found) then
     187     call abort_physic(modname, &
     188                "phyetat0: Failed loading <ZSTD>",1)
     189   endif
     190else
     191  zstd(:)=0.
     192endif ! of if (startphy_file)
     193write(*,*) "phyetat0: <ZSTD> range:", &
     194            minval(zstd), maxval(zstd)
    152195
    153196
    154197! ZSIG
    155 call get_field("ZSIG",zsig,found)
    156 if (.not.found) then
    157   write(*,*) "phyetat0: Failed loading <ZSIG>"
    158   call abort
    159 else
    160   write(*,*) "phyetat0: <ZSIG> range:", &
    161              minval(zsig), maxval(zsig)
    162 endif
     198if (startphy_file) then
     199   call get_field("ZSIG",zsig,found)
     200   if (.not.found) then
     201     call abort_physic(modname, &
     202                "phyetat0: Failed loading <ZSIG>",1)
     203   endif
     204else
     205  zsig(:)=0.
     206endif ! of if (startphy_file)
     207write(*,*) "phyetat0: <ZSIG> range:", &
     208            minval(zsig), maxval(zsig)
    163209
    164210
    165211! ZGAM
    166 call get_field("ZGAM",zgam,found)
    167 if (.not.found) then
    168   write(*,*) "phyetat0: Failed loading <ZGAM>"
    169   call abort
    170 else
    171   write(*,*) "phyetat0: <ZGAM> range:", &
    172              minval(zgam), maxval(zgam)
    173 endif
     212if (startphy_file) then
     213   call get_field("ZGAM",zgam,found)
     214   if (.not.found) then
     215     call abort_physic(modname, &
     216                "phyetat0: Failed loading <ZGAM>",1)
     217   endif
     218else
     219  zgam(:)=0.
     220endif ! of if (startphy_file)
     221write(*,*) "phyetat0: <ZGAM> range:", &
     222            minval(zgam), maxval(zgam)
    174223
    175224
    176225! ZTHE
    177 call get_field("ZTHE",zthe,found)
    178 if (.not.found) then
    179   write(*,*) "phyetat0: Failed loading <ZTHE>"
    180   call abort
    181 else
    182   write(*,*) "phyetat0: <ZTHE> range:", &
     226if (startphy_file) then
     227   call get_field("ZTHE",zthe,found)
     228   if (.not.found) then
     229     call abort_physic(modname, &
     230                "phyetat0: Failed loading <ZTHE>",1)
     231   endif
     232else
     233  zthe(:)=0.
     234endif ! of if (startphy_file)
     235write(*,*) "phyetat0: <ZTHE> range:", &
    183236             minval(zthe), maxval(zthe)
    184 endif
    185 
    186 ! hmons
    187 call get_field("hmons",hmons,found)
    188 if (.not.found) then
    189   write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
    190   if (rdstorm) then
    191     call abort
    192   else
    193     write(*,*) "will continue anyway..."
    194     write(*,*) "because you may not need it."
    195     hmons(:)=0.
    196   end if
    197 else
    198   do ig=1,ngrid
    199     if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
    200   enddo
    201   write(*,*) "phyetat0: <hmons> range:", &
    202              minval(hmons), maxval(hmons)
    203 endif
    204 
    205 ! summit
    206 call get_field("summit",summit,found)
    207 if (.not.found) then
    208   write(*,*) "WARNING: phyetat0: Failed loading <summit>"
    209   if (rdstorm) then
    210     call abort
    211   else
    212     write(*,*) "will continue anyway..."
    213     write(*,*) "because you may not need it."
    214     summit(:)=0.
    215   end if
    216 else
    217   do ig=1,ngrid
    218     if (summit(ig) .eq. -999999.)  summit(ig)=0.
    219   enddo
    220   write(*,*) "phyetat0: <summit> range:", &
    221              minval(summit), maxval(summit)
    222 endif
    223 
    224 ! summit
    225 call get_field("base",base,found)
    226 if (.not.found) then
    227   write(*,*) "WARNING: phyetat0: Failed loading <base>"
    228   if (rdstorm) then
    229     call abort
    230   else
    231     write(*,*) "will continue anyway..."
    232     write(*,*) "because you may not need it."
    233     base(:)=0.
    234   end if
    235 else
    236   do ig=1,ngrid
    237     if (base(ig) .eq. -999999.)  base(ig)=0.
    238   enddo
    239   write(*,*) "phyetat0: <base> range:", &
    240              minval(base), maxval(base)
    241 endif
     237
     238! hmons
     239if (startphy_file) then
     240   call get_field("hmons",hmons,found)
     241   if (.not.found) then
     242     write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
     243     if (rdstorm) then
     244     call abort_physic(modname, &
     245                "phyetat0: Failed loading <hmons>",1)
     246     else
     247       write(*,*) "will continue anyway..."
     248       write(*,*) "because you may not need it."
     249       hmons(:)=0.
     250     end if ! if (rdstorm)
     251   else
     252     do ig=1,ngrid
     253       if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
     254     enddo
     255   endif ! (.not.found)
     256else
     257   hmons(:)=0.
     258endif ! if (startphy_file)
     259write(*,*) "phyetat0: <hmons> range:", &
     260            minval(hmons), maxval(hmons)
     261
     262
     263! summit
     264if (startphy_file) then
     265   call get_field("summit",summit,found)
     266   if (.not.found) then
     267     write(*,*) "WARNING: phyetat0: Failed loading <summit>"
     268     if (rdstorm) then
     269     call abort_physic(modname, &
     270                "phyetat0: Failed loading <summit>",1)
     271     else
     272       write(*,*) "will continue anyway..."
     273       write(*,*) "because you may not need it."
     274       summit(:)=0.
     275     end if
     276   else
     277     do ig=1,ngrid
     278       if (summit(ig) .eq. -999999.)  summit(ig)=0.
     279     enddo
     280   endif ! if (.not.found)
     281else
     282   summit(:)=0. 
     283endif ! if (startphy_file)
     284write(*,*) "phyetat0: <summit> range:", &
     285            minval(summit), maxval(summit)
     286
     287
     288! base
     289if (startphy_file) then
     290   call get_field("base",base,found)
     291   if (.not.found) then
     292     write(*,*) "WARNING: phyetat0: Failed loading <base>"
     293     if (rdstorm) then
     294     call abort_physic(modname, &
     295                "phyetat0: Failed loading <base>",1)
     296     else
     297       write(*,*) "will continue anyway..."
     298       write(*,*) "because you may not need it."
     299       base(:)=0.
     300     end if
     301   else
     302     do ig=1,ngrid
     303       if (base(ig) .eq. -999999.)  base(ig)=0.
     304     enddo
     305   endif ! if(.not.found)
     306else
     307   base(:)=0.
     308endif ! if (startphy_file)
     309write(*,*) "phyetat0: <base> range:", &
     310            minval(base), maxval(base)
     311
    242312 
    243313! Time axis
    244314! obtain timestart from run.def
    245315timestart=-9999 ! default value
    246 call getin("timestart",timestart)
    247 
    248 found=inquire_dimension("Time")
    249 if (.not.found) then
     316call getin_p("timestart",timestart)
     317if (startphy_file) then
     318   found=inquire_dimension("Time")
     319   if (.not.found) then
     320     indextime = 1
     321     write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
     322   else
     323     write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
     324     timelen=inquire_dimension_length("Time")
     325     allocate(time(timelen))
     326     ! load "Time" array:
     327     call get_var("Time",time,found)
     328     if (.not.found) then
     329     call abort_physic(modname, &
     330                "phyetat0: Failed loading <Time>",1)
     331     endif
     332     ! seclect the desired time index
     333     IF (timestart .lt. 0) THEN  ! default: we use the last time value
     334       indextime = timelen
     335     ELSE  ! else we look for the desired value in the time axis
     336       indextime = 0
     337       DO i=1,timelen
     338         IF (abs(time(i) - timestart) .lt. 0.01) THEN
     339           indextime = i
     340           EXIT
     341         ENDIF
     342       ENDDO
     343       IF (indextime .eq. 0) THEN
     344         PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
     345         PRINT*, "Stored times are:"
     346         DO i=1,timelen
     347            PRINT*, time(i)
     348         ENDDO
     349         call abort_physic(modname,"phyetat0: Time error",1)
     350       ENDIF
     351     ENDIF ! of IF (timestart .lt. 0)
     352     ! In startfi the absolute date is day_ini + time0 + time
     353     ! For now on, in the GCM physics, it is day_ini + time0
     354     time0 = time(indextime) + time0
     355     day_ini = day_ini + INT(time0)
     356     time0 = time0 - INT(time0)   
     357     PRINT*, "phyetat0: Selected time ",time(indextime), &
     358             " at index ",indextime
     359     DEALLOCATE(time)
     360   endif ! of if Time not found in file
     361else
    250362  indextime = 1
    251   write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
    252 else
    253   write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
    254   timelen=inquire_dimension_length("Time")
    255   allocate(time(timelen))
    256   ! load "Time" array:
    257   call get_var("Time",time,found)
    258   if (.not.found) then
    259     write(*,*) "phyetat0: Failed loading <Time>"
    260     call abort
    261   endif
    262   ! seclect the desired time index
    263   IF (timestart .lt. 0) THEN  ! default: we use the last time value
    264     indextime = timelen
    265   ELSE  ! else we look for the desired value in the time axis
    266     indextime = 0
    267     DO i=1,timelen
    268       IF (abs(time(i) - timestart) .lt. 0.01) THEN
    269         indextime = i
    270         EXIT
    271       ENDIF
    272     ENDDO
    273     IF (indextime .eq. 0) THEN
    274       PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
    275       PRINT*, "Stored times are:"
    276       DO i=1,timelen
    277          PRINT*, time(i)
    278       ENDDO
    279       CALL abort
    280     ENDIF
    281   ENDIF ! of IF (timestart .lt. 0)
    282   ! In startfi the absolute date is day_ini + time0 + time
    283   ! For now on, in the GCM physics, it is day_ini + time0
    284   time0 = time(indextime) + time0
    285   day_ini = day_ini + INT(time0)
    286   time0 = time0 - INT(time0)
    287        
    288   PRINT*, "phyetat0: Selected time ",time(indextime), &
    289           " at index ",indextime
    290      
    291   DEALLOCATE(time)
    292 endif ! of if Time not found in file
    293 
     363endif ! if (startphy_file)
    294364
    295365! CO2 ice cover
    296 call get_field("co2ice",co2ice,found,indextime)
    297 if (.not.found) then
    298   write(*,*) "phyetat0: Failed loading <co2ice>"
    299   call abort
    300 else
    301   write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
    302              minval(co2ice), maxval(co2ice)
    303 endif
     366if (startphy_file) then
     367   call get_field("co2ice",co2ice,found,indextime)
     368   if (.not.found) then
     369     call abort_physic(modname, &
     370                "phyetat0: Failed loading <co2ice>",1)
     371   endif
     372else
     373   co2ice(:)=0.
     374endif !if (startphy_file)
     375write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
     376            minval(co2ice), maxval(co2ice)
    304377
    305378! Memory of the origin of the co2 particles
    306 call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
    307 if (.not.found) then
    308   write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
    309   mem_Mccn_co2(:,:)=0
    310 else
    311   write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
    312   write(*,*) " <mem_Mccn_co2> range:", &
     379if (startphy_file) then
     380   call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
     381   if (.not.found) then
     382     write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
     383     mem_Mccn_co2(:,:)=0
     384   endif
     385else
     386   mem_Mccn_co2(:,:)=0
     387endif !if (startphy_file)
     388write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
     389write(*,*) " <mem_Mccn_co2> range:", &
    313390             minval(mem_Mccn_co2), maxval(mem_Mccn_co2)
    314 endif
    315 
    316 call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
    317 if (.not.found) then
    318   write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
     391
     392if (startphy_file) then
     393   call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
     394   if (.not.found) then
     395     write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
     396     mem_Nccn_co2(:,:)=0
     397   endif
     398else
    319399  mem_Nccn_co2(:,:)=0
    320 else
    321   write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
    322   write(*,*) " <mem_Nccn_co2> range:", &
     400endif ! if (startphy_file)
     401write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
     402write(*,*) " <mem_Nccn_co2> range:", &
    323403             minval(mem_Nccn_co2), maxval(mem_Nccn_co2)
    324 endif
    325 
    326 call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
    327 if (.not.found) then
    328   write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
    329   mem_Mh2o_co2(:,:)=0
    330 else
    331   write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
    332   write(*,*) " <mem_Mh2o_co2> range:", &
     404
     405if (startphy_file) then
     406   call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
     407   if (.not.found) then
     408     write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
     409     mem_Mh2o_co2(:,:)=0
     410   endif
     411else
     412   mem_Mh2o_co2(:,:)=0
     413endif ! if (startphy_file)
     414write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
     415write(*,*) " <mem_Mh2o_co2> range:", &
    333416             minval(mem_Mh2o_co2), maxval(mem_Mh2o_co2)
    334 endif
    335417
    336418! Dust conversion factor
    337 call get_field("tauscaling",tauscaling,found,indextime)
    338 if (.not.found) then
    339   write(*,*) "phyetat0: <tauscaling> not in file"
    340   tauscaling(:) = 1
    341 else
    342   write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
    343              minval(tauscaling), maxval(tauscaling)
    344 endif
     419if (startphy_file) then
     420   call get_field("tauscaling",tauscaling,found,indextime)
     421   if (.not.found) then
     422     write(*,*) "phyetat0: <tauscaling> not in file"
     423     tauscaling(:) = 1
     424   endif
     425else
     426   tauscaling(:) = 1
     427endif ! if (startphy_file)
     428write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
     429            minval(tauscaling), maxval(tauscaling)
     430
    345431
    346432! dtau: opacity difference between GCM and dust scenario
    347 call get_field("dtau",dtau,found,indextime)
    348 if (.not.found) then
    349   write(*,*) "phyetat0: <dtau> not in file; set to zero"
    350   dtau(:) = 0
    351 else
    352   write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
    353              minval(dtau), maxval(dtau)
    354 endif
     433if (startphy_file) then
     434   call get_field("dtau",dtau,found,indextime)
     435   if (.not.found) then
     436     write(*,*) "phyetat0: <dtau> not in file; set to zero"
     437     dtau(:) = 0
     438   endif
     439else
     440   dtau(:)= 0
     441endif ! if (startphy_file)
     442write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
     443            minval(dtau), maxval(dtau)
    355444
    356445
    357446! Sub-grid cloud fraction
    358 call get_field("totcloudfrac",totcloudfrac,found,indextime)
    359 if (.not.found) then
    360   write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
    361   totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
    362 else
    363   write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
    364              minval(totcloudfrac), maxval(totcloudfrac)
    365 endif
     447if (startphy_file) then
     448   call get_field("totcloudfrac",totcloudfrac,found,indextime)
     449   if (.not.found) then
     450     write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
     451     totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
     452   endif
     453else
     454   totcloudfrac(:)=1.0
     455endif ! if (startphy_file)
     456write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
     457            minval(totcloudfrac), maxval(totcloudfrac)
     458
    366459
    367460! Max vertical velocity in thermals
    368 call get_field("wstar",wstar,found,indextime)
    369 if (.not.found) then
    370   write(*,*) "phyetat0: <wstar> not in file! Set to zero"
    371   wstar(:)=0
    372 else
    373   write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
    374              minval(wstar),maxval(wstar)
    375 endif
     461if (startphy_file) then
     462   call get_field("wstar",wstar,found,indextime)
     463   if (.not.found) then
     464     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
     465     wstar(:)=0
     466   endif
     467else
     468   wstar(:)=0
     469endif ! if (startphy_file)
     470write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
     471            minval(wstar),maxval(wstar)
     472
    376473
    377474! Surface temperature :
    378 call get_field("tsurf",tsurf,found,indextime)
    379 if (.not.found) then
    380   write(*,*) "phyetat0: Failed loading <tsurf>"
    381   call abort
    382 else
    383   write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
    384              minval(tsurf), maxval(tsurf)
    385 endif
     475if (startphy_file) then !tsurf
     476   call get_field("tsurf",tsurf,found,indextime)
     477   if (.not.found) then
     478     call abort_physic(modname, &
     479                "phyetat0: Failed loading <tsurf>",1)
     480   endif
     481else
     482  tsurf(:)=0. ! will be updated afterwards in physiq !
     483endif ! of if (startphy_file)
     484write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
     485            minval(tsurf), maxval(tsurf)
    386486
    387487! Surface albedo
    388 call get_field("albedo",albedo(:,1),found,indextime)
    389 if (.not.found) then
    390   write(*,*) "phyetat0: Failed loading <albedo>"
    391   albedo(:,1)=albedodat(:)
    392 else
    393   write(*,*) "phyetat0: Surface albedo <albedo> range:", &
    394              minval(albedo(:,1)), maxval(albedo(:,1))
    395 endif
     488if (startphy_file) then
     489   call get_field("albedo",albedo(:,1),found,indextime)
     490   if (.not.found) then
     491     write(*,*) "phyetat0: Failed loading <albedo>"
     492     albedo(:,1)=albedodat(:)
     493   endif
     494else
     495   albedo(:,1)=albedodat(:)
     496endif ! of if (startphy_file)
     497write(*,*) "phyetat0: Surface albedo <albedo> range:", &
     498            minval(albedo(:,1)), maxval(albedo(:,1))
    396499albedo(:,2)=albedo(:,1)
    397500
    398501! Surface emissivity
    399 call get_field("emis",emis,found,indextime)
    400 if (.not.found) then
    401   write(*,*) "phyetat0: Failed loading <emis>"
    402   call abort
    403 else
    404   write(*,*) "phyetat0: Surface emissivity <emis> range:", &
    405              minval(emis), maxval(emis)
    406 endif
     502if (startphy_file) then
     503   call get_field("emis",emis,found,indextime)
     504   if (.not.found) then
     505     call abort_physic(modname, &
     506                "phyetat0: Failed loading <emis>",1)
     507   endif
     508else
     509  ! If no startfi file, use parameter surfemis in def file
     510  surfemis=1.0
     511  call getin_p("surfemis",surfemis)
     512  print*,"surfemis",surfemis
     513  emis(:)=surfemis
     514endif ! of if (startphy_file)
     515write(*,*) "phyetat0: Surface emissivity <emis> range:", &
     516            minval(emis), maxval(emis)
    407517
    408518
    409519! surface roughness length (NB: z0 is a common in surfdat_h)
    410 call get_field("z0",z0,found)
    411 if (.not.found) then
    412   write(*,*) "phyetat0: Failed loading <z0>"
    413   write(*,*) 'will use constant value of z0_default:',z0_default
    414   z0(:)=z0_default
    415 else
    416   write(*,*) "phyetat0: Surface roughness <z0> range:", &
    417              minval(z0), maxval(z0)
    418 endif
     520if (startphy_file) then
     521   call get_field("z0",z0,found)
     522   if (.not.found) then
     523     write(*,*) "phyetat0: Failed loading <z0>"
     524     write(*,*) 'will use constant value of z0_default:',z0_default
     525     z0(:)=z0_default
     526   endif
     527else
     528   z0(:)=z0_default
     529endif ! of if (startphy_file)
     530write(*,*) "phyetat0: Surface roughness <z0> range:", &
     531            minval(z0), maxval(z0)
    419532
    420533
    421534! pbl wind variance
    422 call get_field("q2",q2,found,indextime)
    423 if (.not.found) then
    424   write(*,*) "phyetat0: Failed loading <q2>"
    425   call abort
    426 else
    427   write(*,*) "phyetat0: PBL wind variance <q2> range:", &
    428              minval(q2), maxval(q2)
    429 endif
    430 
    431 ! Non-orographic gravity waves
    432 call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
    433 if (.not.found) then
    434    write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
    435    du_nonoro_gwd(:,:)=0.
    436 else
    437    write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
    438    write(*,*) " <du_nonoro_gwd> range:", &
    439           minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
    440 endif
    441 
    442 call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
    443 if (.not.found) then
    444    write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
    445    dv_nonoro_gwd(:,:)=0.
    446 else
    447    write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
    448    write(*,*) " <dv_nonoro_gwd> range:", &
    449           minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
    450 endif
     535if (startphy_file) then
     536   call get_field("q2",q2,found,indextime)
     537   if (.not.found) then
     538     call abort_physic(modname, &
     539                "phyetat0: Failed loading <q2>",1)
     540   endif
     541else
     542  q2(:,:)=0.
     543endif ! of if (startphy_file)
     544write(*,*) "phyetat0: PBL wind variance <q2> range:", &
     545            minval(q2), maxval(q2)
     546
     547! Non-orographic gravity waves
     548if (startphy_file) then
     549   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
     550   if (.not.found) then
     551      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
     552      du_nonoro_gwd(:,:)=0.
     553   endif
     554else
     555du_nonoro_gwd(:,:)=0.
     556endif ! if (startphy_file)
     557write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
     558write(*,*) " <du_nonoro_gwd> range:", &
     559             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
     560
     561if (startphy_file) then
     562   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
     563   if (.not.found) then
     564      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
     565      dv_nonoro_gwd(:,:)=0.
     566   endif
     567else ! ! if (startphy_file)
     568dv_nonoro_gwd(:,:)=0.
     569endif ! if (startphy_file)
     570write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
     571write(*,*) " <dv_nonoro_gwd> range:", &
     572             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
    451573
    452574! tracer on surface
    … …  
    462584      write(*,*) 'iq = ', iq
    463585    endif
    464     call get_field(txt,qsurf(:,iq),found,indextime)
    465     if (.not.found) then
    466       write(*,*) "phyetat0: Failed loading <",trim(txt),">"
    467       write(*,*) "         ",trim(txt)," is set to zero"
     586    if (startphy_file) then
     587       call get_field(txt,qsurf(:,iq),found,indextime)
     588       if (.not.found) then
     589         write(*,*) "phyetat0: Failed loading <",trim(txt),">"
     590         write(*,*) "         ",trim(txt)," is set to zero"
     591       endif
    468592    else
    469       write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
     593      qsurf(:,iq)=0.
     594    endif ! of if (startphy_file)
     595    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
    470596                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
    471     endif
    472   enddo
     597  enddo ! of do iq=1,nq
    473598endif ! of if (nq.ge.1)
    474599
    475 
    476 call get_field("watercap",watercap,found,indextime)
    477 if (.not.found) then
    478   write(*,*) "phyetat0: Failed loading <watercap> : ", &
    479                        "<watercap> is set to zero"
    480   watercap(:)=0
    481 
    482   write(*,*) 'Now transfer negative surface water ice to', &
    483              ' watercap !'
    484   if (nq.ge.1) then
    485    do iq=1,nq
    486     txt=noms(iq)
    487     if (txt.eq."h2o_ice") then
    488       do ig=1,ngrid
    489        if (qsurf(ig,iq).lt.0.0) then
    490           watercap(ig) = qsurf(ig,iq)
    491           qsurf(ig,iq) = 0.0
    492        end if
    493       end do
    494     endif
    495    enddo
    496   endif ! of if (nq.ge.1)
    497 
    498 else
    499   write(*,*) "phyetat0: Surface water ice <watercap> range:", &
    500              minval(watercap), maxval(watercap)
    501 endif
    502 
    503 
    504 ! Call to soil_settings, in order to read soil temperatures,
    505 ! as well as thermal inertia and volumetric heat capacity
    506 call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
    507 
     600if (startphy_file) then
     601   call get_field("watercap",watercap,found,indextime)
     602   if (.not.found) then
     603     write(*,*) "phyetat0: Failed loading <watercap> : ", &
     604                          "<watercap> is set to zero"
     605     watercap(:)=0
     606
     607     write(*,*) 'Now transfer negative surface water ice to', &
     608                ' watercap !'
     609     if (nq.ge.1) then
     610      do iq=1,nq
     611       txt=noms(iq)
     612       if (txt.eq."h2o_ice") then
     613         do ig=1,ngrid
     614          if (qsurf(ig,iq).lt.0.0) then
     615             watercap(ig) = qsurf(ig,iq)
     616             qsurf(ig,iq) = 0.0
     617          end if
     618         end do
     619       endif
     620      enddo
     621     endif ! of if (nq.ge.1)
     622   endif ! of if (.not.found)
     623else
     624   watercap(:)=0
     625endif ! of if (startphy_file)
     626write(*,*) "phyetat0: Surface water ice <watercap> range:", &
     627            minval(watercap), maxval(watercap)
     628 
     629
     630
     631if (startphy_file) then
     632  ! Call to soil_settings, in order to read soil temperatures,
     633  ! as well as thermal inertia and volumetric heat capacity
     634  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
     635endif ! of if (startphy_file)
    508636!
    509637! close file:
    510638!
    511 call close_startphy
     639if (startphy_file) call close_startphy
    512640
    513641end subroutine phyetat0

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2265 r2281  
    270270
    271271!      CHARACTER*80 fichier
    272       INTEGER l,ig,ierr,igout,iq,tapphys
     272      INTEGER l,ig,ierr,igout,iq,tapphys,isoil
    273273
    274274      REAL fluxsurf_lw(ngrid)      !incident LW (IR) surface flux (W.m-2)
    … …  
    479479      REAL hsummit(ngrid)
    480480
     481!      LOGICAL startphy_file
     482
    481483c=======================================================================
    482484
    … …  
    515517     &         mem_Mccn_co2,mem_Nccn_co2,
    516518     &         mem_Mh2o_co2,watercap)
    517 
    518          if (pday.ne.day_ini) then
    519            write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
    520      &                "physics and dynamics"
    521            write(*,*) "dynamics day [pday]: ",pday
    522            write(*,*) "physics day [day_ini]: ",day_ini
    523            call abort_physic("physiq","dynamics day /= physics day",1)
    524          endif
    525 
    526          write (*,*) 'In physiq day_ini =', day_ini
    527519
    528520#else
    … …  
    545537      day_ini = pday
    546538#endif
     539#ifndef MESOSCALE
     540         if (.not.startphy_file) then
     541           ! starting without startfi.nc and with callsoil
     542           ! is not yet possible as soildepth default is not defined
     543           if (callsoil) then
     544              call abort_physic("physiq","callsoil option is not",
     545     &                     "yet available without startfi",1)
     546           endif
     547           ! additionnal "academic" initialization of physics
     548           write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
     549           tsurf(:)=pt(:,1)
     550           write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
     551           do isoil=1,nsoilmx
     552             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
     553           enddo
     554           write(*,*) "Physiq: initializing inertiedat !!"
     555           inertiedat(:,:)=400.
     556           write(*,*) "Physiq: initializing day_ini to pdat !"
     557           day_ini=pday
     558        endif
     559#endif
     560
     561         if (pday.ne.day_ini) then
     562           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
     563     &                "physics and dynamics"
     564           write(*,*) "dynamics day [pday]: ",pday
     565           write(*,*) "physics day [day_ini]: ",day_ini
     566           call abort_physic("physiq","dynamics day /= physics day",1)
     567         endif
     568
     569         write (*,*) 'In physiq day_ini =', day_ini
    547570
    548571c        initialize tracers
    … …  
    654677      ENDIF        !  (end of "if firstcall")
    655678
    656 
    657679c ---------------------------------------------------
    658680c 1.2   Initializations done at every physical timestep:
    … …  
    664686      call update_xios_timestep
    665687#endif     
    666 
    667 c     Initialize various variables
    668688c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    669689      pdv(:,:)=0
    … …  
    705725      zplay(:,:) = pplay(:,:)
    706726      ps(:) = pplev(:,1)
    707 
    708727
    709728c     Compute geopotential at interlayers
    … …  
    739758         ENDDO
    740759      ENDDO
    741 
    742760#ifndef MESOSCALE
    743761c-----------------------------------------------------------------------
    … …  
    768786c    2. Compute radiative tendencies :
    769787c------------------------------------
    770 
    771788
    772789      IF (callrad) THEN
    … …  
    931948     &           zls*180./pi,latitude(igout)*180/pi,
    932949     &                       longitude(igout)*180/pi
     950
    933951           write(*,'(" tauref(",f4.0," Pa) =",f9.6,
    934952     &             " tau(",f4.0," Pa) =",f9.6)')
    … …  
    11611179c    -------------------------------------------
    11621180           if (dustinjection.gt.0) then
     1181
    11631182             CALL compute_dtau(ngrid,nlayer,
    11641183     &                          zday,pplev,tauref,
    … …  
    12291248             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
    12301249         ENDIF
    1231 
    12321250c        Calling vdif (Martian version WITH CO2 condensation)
    12331251         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
    … …  
    12561274            ENDDO
    12571275          ENDDO
    1258 
    12591276          if (tracer) then
    12601277           DO iq=1, nq
    … …  
    35403557
    35413558      icount=icount+1
    3542        
    35433559
    35443560      END SUBROUTINE physiq

  • trunk/LMDZ.MARS/libf/phymars/tabfi.F

    r1543 r2281  
    4343c
    4444c=======================================================================
     45! to use  'getin_p'
     46      use ioipsl_getin_p_mod, only: getin_p
    4547
    4648      use comsoil_h, only: volcapa ! soil volumetric heat capacity
    … …  
    8183      CHARACTER modif*20
    8284      LOGICAL :: found
     85      CHARACTER(len=5) :: modname="tabfi"
    8386
    8487      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
    … …  
    8891c-----------------------------------------------------------------------
    8992      IF (nid.eq.0) then
     93
     94      ! to avoid further issues with writing
     95      tab_cntrl(:)=0
     96      lmax=0
     97      day_ini=0
     98      time = 0
    9099 
    91100c Reference pressure
    … …  
    150159       call get_var("controle",tab_cntrl,found)
    151160       if (.not.found) then
    152          write(*,*)"tabfi: Failed reading <controle> array"
    153          call abort
     161         call abort_physic(modname,
     162     &        "tabfi: Failed reading <controle> array",1)
    154163       else
    155164         write(*,*)'tabfi: tab_cntrl',tab_cntrl
Note: See TracChangeset for help on using the changeset viewer.