Ignore:
Timestamp:
Sep 26, 2023, 3:57:38 PM (15 months ago)
Author:
jbclement
Message:

Mars PEM:
Minor changes concerning the form of the code in the PEM.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
8 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90

    r3039 r3050  
    6363     present_surf.GT.ini_surf*(1+water_ice_criterion)) then
    6464    STOPPING=.TRUE.
    65     print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
    66     print *, "Current surface of water ice sublimating=", present_surf
    67     print *, "Initial surface of water ice sublimating=", ini_surf
    68     print *, "Percentage of change accepted=", water_ice_criterion*100
    69     print *, "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion))
     65    write(*,*) "Reason of stopping : The surface of water ice sublimating reach the threshold:"
     66    write(*,*) "Current surface of water ice sublimating=", present_surf
     67    write(*,*) "Initial surface of water ice sublimating=", ini_surf
     68    write(*,*) "Percentage of change accepted=", water_ice_criterion*100
     69    write(*,*) "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion))
    7070  endif
    7171
    … …  
    129129     present_surf.GT.ini_surf*(1+co2_ice_criterion)) then
    130130    STOPPING_ice=.TRUE.
    131     print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
    132     print *, "Current surface of co2 ice sublimating=", present_surf
    133     print *, "Initial surface of co2 ice sublimating=", ini_surf
    134     print *, "Percentage of change accepted=", co2_ice_criterion*100
    135     print *, "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion))
     131    write(*,*) "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
     132    write(*,*) "Current surface of co2 ice sublimating=", present_surf
     133    write(*,*) "Initial surface of co2 ice sublimating=", ini_surf
     134    write(*,*) "Percentage of change accepted=", co2_ice_criterion*100
     135    write(*,*) "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion))
    136136  endif
    137137
    … …  
    143143     global_ave_press_new.GT.global_ave_press_GCM*(1+ps_criterion)) then
    144144    STOPPING_ps=.TRUE.
    145     print *, "Reason of stopping : The global pressure reach the threshold:"
    146     print *, "Current global pressure=", global_ave_press_new
    147     print *, "GCM global pressure=", global_ave_press_GCM
    148     print *, "Percentage of change accepted=", ps_criterion*100
    149     print *, "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion))
     145    write(*,*) "Reason of stopping : The global pressure reach the threshold:"
     146    write(*,*) "Current global pressure=", global_ave_press_new
     147    write(*,*) "GCM global pressure=", global_ave_press_GCM
     148    write(*,*) "Percentage of change accepted=", ps_criterion*100
     149    write(*,*) "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion))
    150150  endif
    151151

  • trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90

    r3039 r3050  
    9494     enddo
    9595   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then
    96     print *, "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
    97     print *, "Tendencies on ice sublimating=", neg_tend
    98     print *, "Tendencies on ice increasing=", pos_tend
    99     print *, "This can be due to the absence of water ice in the PCM run!!"
     96    write(*,*) "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
     97    write(*,*) "Tendencies on ice sublimating=", neg_tend
     98    write(*,*) "Tendencies on ice increasing=", pos_tend
     99    write(*,*) "This can be due to the absence of water ice in the PCM run!!"
    100100    call criterion_waterice_stop(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,qsurf(:,:)*0.)
    101101    do i=1,ngrid

  • trunk/LMDZ.COMMON/libf/evolution/init_pem1d_mod.F90

    r3049 r3050  
    1 module init_phys_1d_mod
     1module init_pem1d_mod
    22
    33implicit none
    … …  
    77#ifdef CPP_1D
    88
    9 subroutine init_phys_1d(llm_in,nqtot_in,v,u,temp,q,psurf,time,ap_out,bp_out)
    10 
    11 !-------------- init_phys_1d -------------
     9subroutine init_pem1d(llm_in,nqtot_in,u,v,temp,q,psurf,time,ap_out,bp_out)
     10
     11!-------------- init_pem1d -------------
    1212!
    1313! This function is a copy of the test_phys_1d.F90 initialisation part.
    … …  
    2121      USE ioipsl_getincom, only: getin
    2222      use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
    23       use time_phylmdz_mod, only: daysec, dtphys, day_step, &
    24                                  ecritphy, iphysiq
    25       use planete_h, only: year_day, periheli, aphelie, peri_day,&
    26                           obliquit, emin_turb, lmixmin
     23      use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy, iphysiq
     24      use planete_h, only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin
    2725      use surfdat_h, only: albedodat, z0_default, emissiv, emisice,&
    2826                          albedice, iceradius, dtemisice, z0,&
    2927                          zmea, zstd, zsig, zgam, zthe, phisfi,&
    3028                          watercaptag, watercap, hmons, summit, base,&
    31                            tsurf, emis,qsurf       
     29                          tsurf, emis,qsurf       
    3230      use infotrac, only: nqtot, tname, nqperes,nqfils
    3331      use regular_lonlat_mod, only: init_regular_lonlat
    … …  
    6967!-------------- OUTPUT VARIABLES -------------
    7068
    71       REAL,INTENT(OUT):: u(nlayer),v(nlayer) ! zonal, meridional wind
    72       REAL,INTENT(OUT):: temp(nlayer)   ! temperature at the middle of the layers
    73       REAL,INTENT(OUT) :: q(2,llm_in,nqtot_in) ! tracer mixing ratio (e.g. kg/kg)
    74       REAL,INTENT(OUT):: psurf(2)
    75       REAL,INTENT(OUT):: time             ! time (0<time<1 ; time=0.5 a midi)
    76       REAL,INTENT(OUT) :: ap_out(llm_in+1),bp_out(llm_in+1)
     69      real, intent(out) :: u(nlayer), v(nlayer) ! zonal, meridional wind
     70      real, intent(out) :: temp(nlayer)         ! temperature at the middle of the layers
     71      real, intent(out) :: q(2,llm_in,nqtot_in) ! tracer mixing ratio (e.g. kg/kg)
     72      real, intent(out) :: psurf(2)
     73      real, intent(out) :: time                 ! time (0<time<1 ; time=0.5 a midi)
     74      real, intent(out) :: ap_out(llm_in + 1), bp_out(llm_in + 1)
    7775
    7876
    … …  
    10098      REAL tmp1(0:nlayer),tmp2(0:nlayer)
    10199      INTEGER flagthermo,flagh2o
    102       REAL atm_wat_profile
    103       REAL,ALLOCATABLE :: zqsat(:) ! useful for (atm_wat_profile=2)
     100      real atm_wat_profile, atm_wat_tau
     101      real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
    104102
    105103      INTEGER :: ierr,iq,ilayer,ilevel,isoil
    … …  
    112110      real :: flux_geo_tmp
    113111
     112!   RV & JBC: Use of "start1D.txt" and "startfi.nc" files
     113      logical                :: therestart1D, therestartfi
     114      character(len = 30)    :: header
     115
    114116      real,allocatable :: qdyn(:,:,:,:),psdyn(:,:)
    115117
    116118
    117 ! check if 'run.def' file is around (otherwise reading parameters
    118 ! from callphys.def via getin() routine won't work.
    119       open(99,file='run.def',status='old',form='formatted',&
    120           iostat=ierr)
    121       if (ierr.ne.0) then
    122         write(*,*) 'Cannot find required file "run.def"'
    123         write(*,*) '  (which should contain some input parameters'
    124         write(*,*) '   along with the following line:'
    125         write(*,*) '   INCLUDEDEF=callphys.def'
    126         write(*,*) '   )'
    127         write(*,*) ' ... might as well stop here ...'
    128         stop
    129       else
    130         close(99)
    131       endif
     119! Checking if the file "start1D.txt" exists
     120inquire(file ='start1D.txt',exist = therestart1D)
     121if (.not. therestart1D) then
     122    write(*,*) 'There is no "start1D.txt" file!'
     123    write(*,*) 'It is required to initialize the 1D PEM.'
     124    stop
     125endif
     126inquire(file ='startfi.nc',exist = therestartfi)
     127if (.not. therestartfi) then
     128    write(*,*) 'There is no "startfi.nc" file!'
     129    write(*,*) 'Initialization is done with default values.'
     130endif
    132131
    133132! ------------------------------------------------------
    134133!  Prescribed constants to be set here
    135134! ------------------------------------------------------
    136 
    137 !      if(.not. startfile_1D ) then
    138 
    139135      pi=2.E+0*asin(1.E+0)
    140136
    … …  
    175171      dtemisice(2) = 2.          ! time scale for snow metamorphism (south)
    176172
    177 !      endif ! .not. startfile_1D
    178 
    179173!     mesh surface (not a very usefull quantity in 1D)
    180174!     ----------------------------------------------------
    181175      cell_area(1)=1.E+0
    182176
    183 
    184 ! check if there is a 'traceur.def' file
    185 ! and process it.
    186       ! load tracer names from file 'traceur.def'
    187         open(90,file='traceur.def',status='old',form='formatted',&
    188             iostat=ierr)
    189         if (ierr.ne.0) then
    190           write(*,*) 'Cannot find required file "traceur.def"'
    191           write(*,*) ' If you want to run with tracers, I need it'
    192           write(*,*) ' ... might as well stop here ...'
    193           stop
    194         else
    195           nqtot=nqtot_in ! set value of nqtot (in infotrac module) as nq
    196           nq=nqtot_in
     177! check if there is a 'traceur.def' file and process it
     178! load tracer names from file 'traceur.def'
     179    open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
     180    if (ierr /= 0) then
     181        write(*,*) 'Cannot find required file "traceur.def"'
     182        write(*,*) ' If you want to run with tracers, I need it'
     183        write(*,*) ' ... might as well stop here ...'
     184        stop
     185    else
     186        write(*,*) "pem1d: Reading file traceur.def"
     187        ! read number of tracers:
     188        read(90,*,iostat = ierr) nq
     189        write(*,*) "nq",nq
     190        nqtot = nq ! set value of nqtot (in infotrac module) as nq
     191        if (ierr /= 0) then
     192            write(*,*) "pem1d: error reading number of tracers"
     193            write(*,*) "   (first line of traceur.def) "
     194            stop
    197195        endif
    198         ! allocate arrays:
    199         allocate(tname(nq))
    200         allocate(qsurf(1,1,nq))
    201         allocate(tnom_transp(nq))
    202        
    203         ! read tracer names from file traceur.def
    204         do iq=1,nq
    205           read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
    206           if (ierr.ne.0) then
    207             write(*,*) 'testphys1d: error reading tracer names...'
     196        if (nq < 1) then
     197            write(*,*) "pem1d: error number of tracers"
     198            write(*,*) "is nq=",nq," but must be >=1!"
    208199            stop
    209           endif
    210           ! if format is tnom_0, tnom_transp (isotopes)
    211           read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq)
    212           if (ierr0.ne.0) then
     200        endif
     201    endif
     202    ! allocate arrays:
     203    allocate(tname(nq))
     204    allocate(qsurf(1,1,nq))
     205    allocate(tnom_transp(nq))
     206
     207    ! read tracer names from file traceur.def
     208    do iq = 1,nq
     209        read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
     210        if (ierr /= 0) then
     211            write(*,*) 'init_pem1d: error reading tracer names...'
     212            stop
     213        endif
     214        ! if format is tnom_0, tnom_transp (isotopes)
     215        read(line,*,iostat = ierr0) tname(iq),tnom_transp(iq)
     216        if (ierr0 /= 0) then
    213217            read(line,*) tname(iq)
    214218            tnom_transp(iq)='air'
    215           endif
    216 
    217         enddo
    218         close(90)
     219        endif
     220    enddo
     221    close(90)
    219222
    220223       ! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
    … …  
    225228       if (tnom_transp(iq) == 'air') then
    226229         ! ceci est un traceur père
    227          WRITE(*,*) 'Le traceur',iq,', appele ',&
    228                trim(tname(iq)),', est un pere'
     230         WRITE(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un pere'
    229231         nqperes=nqperes+1
    230232       else !if (tnom_transp(iq) == 'air') then
    231233         ! ceci est un fils. Qui est son père?
    232          WRITE(*,*) 'Le traceur',iq,', appele ',&
    233                      trim(tname(iq)),', est un fils'
     234         WRITE(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un fils'
    234235         continu=.true.
    235236         ipere=1
    … …  
    237238           if (tnom_transp(iq) .eq. tname(ipere)) then
    238239             ! Son père est ipere
    239              WRITE(*,*) 'Le traceur',iq,'appele ',&
    240         trim(tname(iq)),' est le fils de ',&
    241         ipere,'appele ',trim(tname(ipere))
     240             WRITE(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),' est le fils de ',ipere,'appele ',trim(tname(ipere))
    242241             nqfils(ipere)=nqfils(ipere)+1         
    243242             continu=.false.
    … …  
    245244             ipere=ipere+1
    246245             if (ipere.gt.nqtot) then
    247                  WRITE(*,*) 'Le traceur',iq,'appele ',&
    248                 trim(tname(iq)),', est orpelin.'
    249                  CALL abort_gcm('infotrac_init',&
    250                        'Un traceur est orphelin',1)
     246                 WRITE(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.'
     247                 CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
    251248             endif !if (ipere.gt.nqtot) then
    252249           endif !if (tnom_transp(iq) == tnom_0(ipere)) then
    … …  
    257254       WRITE(*,*) 'nqfils=',nqfils
    258255
    259         ! initialize tracers here:
    260         write(*,*) "testphys1d: initializing tracers"
    261         call read_profile(nq, nlayer, qsurf, q(1,:,:))
    262         q(2,:,:)=q(1,:,:)
     256        ! Initialize tracers here:
     257        write(*,*) "init_pem1d: initializing tracers"
     258        do iq = 1,nq
     259            open(3,file = 'start1D.txt',status = "old",action = "read")
     260            read(3,*) header, qsurf(1,1,iq), (q(1,ilayer,iq), ilayer = 1,nlayer)
     261            if (tname(iq) /= trim(header)) then
     262                write(*,*) 'Tracer names not compatible for initialization with "start1D.txt"!'
     263                stop
     264            endif
     265        enddo
     266        q(2,:,:) = q(1,:,:)
    263267
    264268#ifdef CPP_XIOS
    265       call init_physics_distribution(regular_lonlat,4,&
    266                                     1,1,1,nlayer,COMM_LMDZ)
     269    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ)
    267270#else
    268       call init_physics_distribution(regular_lonlat,4,&
    269                                     1,1,1,nlayer,1)
     271    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1)
    270272#endif
    271 
    272       call open_startphy("startfi_evol.nc")
    273       call get_var("controle",tab_cntrl,found)
    274        if (.not.found) then
    275          call abort_physic("open_startphy", &
    276              "tabfi: Failed reading <controle> array",1)
    277        else
    278          write(*,*)'tabfi: tab_cntrl',tab_cntrl
    279        endif
    280        day0 = tab_cntrl(3)
    281        day=float(day0)
    282 
    283        call get_var("Time",time,found)
    284 
    285        call close_startphy
    286 
    287273
    288274!  Discretization (Definition of grid and time steps)
    289275!  --------------
    290 !
    291276      nlevel=nlayer+1
    292277      nsoil=nsoilmx
    293278
    294279      day_step=48 ! default value for day_step
    295       PRINT *,'Number of time steps per sol ?'
     280      write(*,*)'Number of time steps per sol ?'
    296281      call getin("day_step",day_step)
    297282      write(*,*) " day_step = ",day_step
    … …  
    300285
    301286      ndt=10 ! default value for ndt
    302       PRINT *,'Number of sols to run ?'
     287      write(*,*)'Number of sols to run ?'
    303288      call getin("ndt",ndt)
    304289      write(*,*) " ndt = ",ndt
    … …  
    306291      dayn=day0+ndt
    307292      ndt=ndt*day_step     
    308       dtphys=daysec/day_step 
     293      dtphys=daysec/day_step
    309294
    310295! Imposed surface pressure
    311296! ------------------------------------
    312 !
    313 
    314       psurf(1)=610. ! default value for psurf
    315       PRINT *,'Surface pressure (Pa) ?'
    316       call getin("psurf",psurf(1))
     297      psurf(1) = 610. ! Default value for psurf
     298      write(*,*) 'Surface pressure (Pa)?'
     299      if (.not. therestart1D) then
     300          call getin("psurf",psurf)
     301      else
     302          read(3,*) header, psurf(1)
     303      endif
    317304      write(*,*) " psurf = ",psurf(1)
    318       psurf(2)=psurf(1)
     305      psurf(2) = psurf(1)
    319306
    320307! Reference pressures
    … …  
    329316      write(*,*) " tauvis = ",tauvis
    330317
    331 ! Orbital parameters
    332 ! ------------------
    333 
    334318!  latitude/longitude
    335319!  ------------------
    336320      latitude(1)=0 ! default value for latitude
    337       PRINT *,'latitude (in degrees) ?'
     321      write(*,*)'latitude (in degrees) ?'
    338322      call getin("latitude",latitude(1))
    339323      write(*,*) " latitude = ",latitude
    … …  
    344328!  some initializations (some of which have already been
    345329!  done above!) and loads parameters set in callphys.def
    346 !  and allocates some arrays 
    347 !Mars possible matter with dtphys in input and include!!!
     330!  and allocates some arrays
     331! Mars possible matter with dtphys in input and include!!!
    348332! Initializations below should mimick what is done in iniphysiq for 3D GCM
    349333      call init_interface_dyn_phys
    … …  
    368352      ! ovverride iphysiq value that has been set by conf_phys
    369353      if (iphysiq/=1) then
    370         write(*,*) "testphys1d: setting iphysiq=1"
     354        write(*,*) "init_pem1d: setting iphysiq=1"
    371355        iphysiq=1
    372       endif
    373 
    374 !  Initialize albedo / soil thermal inertia
    375 !  ----------------------------------------
    376 !
    377 
     356      endif
    378357
    379358! Initialize local slope parameters (only matters if "callslope"
    … …  
    390369      def_slope(1)=-90 ! minimum slope angle
    391370      def_slope(2)=90 ! maximum slope angle
    392       subslope_dist(1,1)=1 ! fraction of subslopes in mesh 
     371      subslope_dist(1,1)=1 ! fraction of subslopes in mesh
    393372!
    394373!  for the gravity wave scheme
    395374!  ---------------------------------
    396 !
    397375      zmea(1)=0.E+0
    398376      zstd(1)=0.E+0
    … …  
    405383! 
    406384      hmons(1)=0.E+0
    407       PRINT *,'hmons is initialized to ',hmons(1)
     385      write(*,*)'hmons is initialized to ',hmons(1)
    408386      summit(1)=0.E+0
    409       PRINT *,'summit is initialized to ',summit(1)
     387      write(*,*)'summit is initialized to ',summit(1)
    410388      base(1)=0.E+0
    411389!
    412390!  Default values initializing the coefficients calculated later
    413391!  ---------------------------------
    414 ! 
    415392      tauscaling(1)=1. ! calculated in aeropacity_mod.F
    416393      totcloudfrac(1)=1. ! calculated in watercloud_mod.F     
    … …  
    428405!    geostrophic wind
    429406      gru=10. ! default value for gru
    430       PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
     407      write(*,*)'zonal eastward component of the geostrophic wind (m/s) ?'
    431408      call getin("u",gru)
    432409      write(*,*) " u = ",gru
    433410      grv=0. !default value for grv
    434       PRINT *,'meridional northward component of the geostrophic',&
     411      write(*,*)'meridional northward component of the geostrophic',&
    435412     ' wind (m/s) ?'
    436413      call getin("v",grv)
    437414      write(*,*) " v = ",grv
    438415
    439 !     Initialize winds  for first time step
    440       DO ilayer=1,nlayer
    441          u(ilayer)=gru
    442          v(ilayer)=grv
    443          w(ilayer)=0 ! default: no vertical wind
    444       ENDDO
     416!     Initialize winds for first time step
     417      read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
     418      read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
     419      w = 0. ! default: no vertical wind
    445420
    446421!     Initialize turbulente kinetic energy
    … …  
    459434      enddo
    460435      if (igcm_co2==0) then
    461         write(*,*) "testphys1d error, missing co2 tracer!"
     436        write(*,*) "init_pem1d error, missing co2 tracer!"
    462437        stop
    463438      endif
    … …  
    465440!  Compute pressures and altitudes of atmospheric levels
    466441!  ----------------------------------------------------------------
    467 
    468442!    Vertical Coordinates
    469443!    """"""""""""""""""""
    470444      hybrid=.true.
    471       PRINT *,'Hybrid coordinates ?'
     445      write(*,*)'Hybrid coordinates ?'
    472446      call getin("hybrid",hybrid)
    473447      write(*,*) " hybrid = ", hybrid
    … …  
    487461
    488462      DO ilayer=1,nlayer
    489         zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))&
    490         /g
     463        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))/g
    491464      ENDDO
    492465
    … …  
    504477      call profile(nlayer+1,tmp1,tmp2)
    505478
    506       tsurf=tmp2(0)
    507       DO ilayer=1,nlayer
    508         temp(ilayer)=tmp2(ilayer)
    509       ENDDO
    510      
     479      read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer)
     480      close(3)
    511481
    512482! Initialize soil properties and temperature
    … …  
    527497      enddo
    528498
    529 !    Initialize traceurs
    530 !    ---------------------------
    531 
     499! Initialize traceurs
     500! ---------------------------
    532501      if (photochem.or.callthermos) then
    533502         write(*,*) 'Initializing chemical species'
    … …  
    548517         qdyn(1,1,1:llm,1:nq)=q(1,1:llm,1:nq)
    549518         psdyn(1:2,1)=psurf(1)
    550          call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,&
    551                               flagh2o,flagthermo)
     519         call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,flagh2o,flagthermo)
    552520         q(1,1:llm,1:nq)=qdyn(1,1,1:llm,1:nq)
    553521      endif
    … …  
    556524! --------------------------------------------------
    557525      watercaptag(1)=.false. ! Default: no water ice reservoir
    558       print *,'Water ice cap on ground ?'
     526      write(*,*)'Water ice cap on ground ?'
    559527      call getin("watercaptag",watercaptag)
    560528      write(*,*) " watercaptag = ",watercaptag
    … …  
    563531! ---------------------------------------------------
    564532      ! Adding an option to force atmospheric water values JN
    565       atm_wat_profile=-1 ! Default: free atm wat profile
    566       print *,'Force atmospheric water vapor profile ?'
    567       call getin("atm_wat_profile",atm_wat_profile)
    568       write(*,*) "atm_wat_profile = ", atm_wat_profile
    569       if (atm_wat_profile.EQ.-1) then
    570         write(*,*) "Free atmospheric water vapor profile"
    571         write(*,*) "Total water is conserved on the column"
    572       else if (atm_wat_profile.EQ.0) then
    573         write(*,*) "Dry atmospheric water vapor profile"
    574       else if ((atm_wat_profile.GT.0).and.(atm_wat_profile.LE.1)) then
    575         write(*,*) "Prescribed atmospheric water vapor MMR profile"
    576         write(*,*) "Unless it reaches saturation (maximal value)"
    577         write(*,*) "MMR chosen : ", atm_wat_profile
    578       endif
    579 
    580       ap_out=ap
    581       bp_out=bp
    582 
    583 end subroutine init_phys_1d
     533      atm_wat_profile = -1. ! Default: free atm wat profile
     534      if (water) then
     535          write(*,*)'Force atmospheric water vapor profile?'
     536          call getin('atm_wat_profile',atm_wat_profile)
     537          write(*,*) 'atm_wat_profile = ', atm_wat_profile
     538          if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
     539              write(*,*) 'Free atmospheric water vapor profile'
     540              write(*,*) 'Total water is conserved in the column'
     541          else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
     542              write(*,*) 'Dry atmospheric water vapor profile'
     543          else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
     544              write(*,*) 'Prescribed atmospheric water vapor profile'
     545              write(*,*) 'Unless it reaches saturation (maximal value)'
     546          else
     547              write(*,*) 'Water vapor profile value not correct!'
     548              stop
     549          endif
     550      endif
     551
     552! Check if the atmospheric water profile relaxation is specified
     553! --------------------------------------------------------------
     554      ! Adding an option to relax atmospheric water values JBC
     555      atm_wat_tau = -1. ! Default: no time relaxation
     556      if (water) then
     557          write(*,*) 'Relax atmospheric water vapor profile?'
     558          call getin('atm_wat_tau',atm_wat_tau)
     559          write(*,*) 'atm_wat_tau = ', atm_wat_tau
     560          if (atm_wat_tau < 0.) then
     561              write(*,*) 'Atmospheric water vapor profile is not relaxed.'
     562          else
     563              if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
     564                  write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile
     565                  write(*,*) 'Unless it reaches saturation (maximal value)'
     566              else
     567                  write(*,*) 'Reference atmospheric water vapor profile not known!'
     568                  write(*,*) 'Please, specify atm_wat_profile'
     569                  stop
     570              endif
     571          endif
     572      endif
     573
     574      ap_out = ap
     575      bp_out = bp
     576
     577end subroutine init_pem1d
    584578
    585579#endif
    586580
    587 end module init_phys_1d_mod
     581end module init_pem1d_mod

  • trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90

    r2855 r3050  
    316316    ELSE
    317317      IF (.NOT. tmp_found) THEN
    318         PRINT*, 'get_field_rgen: Field <'//field_name//'> not found'
     318        write(*,*) 'get_field_rgen: Field <'//field_name//'> not found'
    319319        CALL abort_physic("get_field_rgen","Failed to get field",1)
    320320      ENDIF
    … …  
    329329         IF (ierr/=NF90_NOERR) THEN
    330330           ! La variable exist dans le fichier mais la lecture a echouee.
    331            PRINT*, 'get_field_rgen: Failed reading <'//field_name//'>'
     331           write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'
    332332
    333333!           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
    … …  
    336336!              ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo))
    337337!              IF (ierr/=NF90_NOERR) THEN
    338 !                 PRINT*, 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
     338!                 write(*,*) 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
    339339!                 CALL abort
    340340!              ELSE
    341 !                 PRINT*, 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
     341!                 write(*,*) 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
    342342!              END IF
    343343!           ELSE
    … …  
    436436        ierr=NF90_GET_VAR(nid_start,varid,var)
    437437        IF (ierr/=NF90_NOERR) THEN
    438           PRINT*, 'phyetat0: Failed loading <'//trim(var_name)//'>'
     438          write(*,*) 'phyetat0: Failed loading <'//trim(var_name)//'>'
    439439          CALL abort_physic("get_var_rgen","Failed to read variable",1)
    440440        ENDIF
    … …  
    456456    ELSE
    457457      IF (.NOT. tmp_found) THEN
    458         PRINT*, 'phyetat0: Variable <'//trim(var_name)//'> not found'
     458        write(*,*) 'phyetat0: Variable <'//trim(var_name)//'> not found'
    459459        CALL abort_physic("get_var_rgen","Failed to read variable",1)
    460460      ENDIF
    … …  
    900900
    901901      ELSE
    902         PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
     902        write(*,*) "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
    903903        write(*,*) "  field_size =",field_size
    904904        CALL abort_physic("put_field_rgen","wrong field dimension",1)
    … …  
    10261026        idim1d=idim9
    10271027      ELSE
    1028         PRINT *, "put_var_rgen error : wrong dimension"
     1028        write(*,*) "put_var_rgen error : wrong dimension"
    10291029        write(*,*) "  var_size =",var_size
    10301030        CALL abort_physic("get_var_rgen","Wrong variable dimension",1)

  • trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90

    r2980 r3050  
    6666  CALL err(NF90_CLOSE(fID),"close",fichnom)
    6767
    68   print *, "The number of timestep of the PCM run data=", timelen
     68  write(*,*) "The number of timestep of the PCM run data=", timelen
    6969
    7070  CONTAINS

  • trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90

    r3047 r3050  
    7979close(73)
    8080if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then
    81     write(*,*) 'The current year could not be find in the file "obl_ecc_lsp.asc".'
     81    write(*,*) 'The current year could not be found in the file "obl_ecc_lsp.asc".'
    8282    stop
    8383else
    … …  
    196196enddo
    197197
    198 if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be find in the file "obl_ecc_lsp.asc".'
    199 if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be find in the file "obl_ecc_lsp.asc".'
    200 if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be find in the file "obl_ecc_lsp.asc".'
     198if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".'
     199if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".'
     200if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".'
    201201
    202202write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3042 r3050  
    9696    use physics_distribution_mod, only: init_physics_distribution
    9797    use mod_grid_phy_lmdz,        only: regular_lonlat
    98     use init_phys_1d_mod,         only: init_phys_1d
     98    use init_pem1d_mod,           only: init_pem1d
    9999#endif
    100100
    … …  
    316316    allocate(ap(nlayer + 1))
    317317    allocate(bp(nlayer + 1))
    318     call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp)
    319     pi = 2.*asin(1.) 
    320     g = 3.72 
     318    call init_pem1d(llm,nqtot,ucov,vcov,teta,q,ps,time_0,ap,bp)
     319    pi = 2.*asin(1.)
     320    g = 3.72
    321321    nsplit_phys = 1
    322322#endif
    … …  
    663663!------------------------
    664664#ifndef CPP_STD
    665      call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
     665    call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    666666#else
    667      call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
     667    call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    668668#endif
    669669
    … …  
    691691        enddo
    692692    enddo
    693     write(*,*) 'Global average pressure old time step',global_ave_press_old
    694693    call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
    695694     
    … …  
    697696        do i = 1,ngrid
    698697            global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    699         enddo 
    700         write(*,*) 'Global average pressure old time step',global_ave_press_old
    701         write(*,*) 'Global average pressure new time step',global_ave_press_new
    702     endif
     698        enddo
     699    endif
     700    write(*,*) 'Global average pressure old time step',global_ave_press_old
     701    write(*,*) 'Global average pressure new time step',global_ave_press_new
    703702
    704703! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
    … …  
    11801179call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
    11811180             float(day_ini),0.,nslope,def_slope,subslope_dist)
    1182 
    1183 
    11841181call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
    11851182             TI_PEM, porefillingice_depth,porefillingice_thickness,         &

  • trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90

    r2985 r3050  
    8080      B=1/m_noco2
    8181
    82   print *, "Opening ", fichnom, "..."
     82  write(*,*) "Opening ", fichnom, "..."
    8383
    8484!  Open initial state NetCDF file
    … …  
    8686  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    8787
    88      print *, "Downloading data for vmr co2..."
     88     write(*,*) "Downloading data for vmr co2..."
    8989
    9090  CALL get_var3("co2_layer1"   ,q_co2_dyn)
    9191
    92      print *, "Downloading data for vmr co2 done"
    93      print *, "Downloading data for vmr h20..."
     92     write(*,*) "Downloading data for vmr co2 done"
     93     write(*,*) "Downloading data for vmr h20..."
    9494
    9595  CALL get_var3("h2o_layer1"   ,q_h2o_dyn)
    9696
    97      print *, "Downloading data for vmr h2o done"
    98      print *, "Downloading data for surface pressure ..."
     97     write(*,*) "Downloading data for vmr h2o done"
     98     write(*,*) "Downloading data for surface pressure ..."
    9999
    100100  CALL get_var3("ps"   ,ps_GCM)
    101101
    102      print *, "Downloading data for surface pressure done"
    103      print *, "nslope=", nslope
     102     write(*,*) "Downloading data for surface pressure done"
     103     write(*,*) "nslope=", nslope
    104104
    105105if(nslope.gt.1) then
    106106
    107      print *, "Downloading data for co2ice_slope ..."
     107     write(*,*) "Downloading data for co2ice_slope ..."
    108108
    109109DO islope=1,nslope
    … …  
    112112ENDDO
    113113
    114      print *, "Downloading data for co2ice_slope done"
    115      print *, "Downloading data for h2o_ice_s_slope ..."
     114     write(*,*) "Downloading data for co2ice_slope done"
     115     write(*,*) "Downloading data for h2o_ice_s_slope ..."
    116116
    117117DO islope=1,nslope
    … …  
    120120ENDDO
    121121
    122      print *, "Downloading data for h2o_ice_s_slope done"
     122     write(*,*) "Downloading data for h2o_ice_s_slope done"
    123123
    124124#ifndef CPP_STD
    125125
    126      print *, "Downloading data for watercap_slope ..."
     126     write(*,*) "Downloading data for watercap_slope ..."
    127127DO islope=1,nslope
    128128       write(num,fmt='(i2.2)') islope
    … …  
    130130!        watercap_slope(:,:,:,:)= 0.
    131131ENDDO           
    132      print *, "Downloading data for watercap_slope done"
     132     write(*,*) "Downloading data for watercap_slope done"
    133133
    134134#endif
    135135   
    136  print *, "Downloading data for tsurf_slope ..."
     136 write(*,*) "Downloading data for tsurf_slope ..."
    137137
    138138DO islope=1,nslope
    … …  
    141141ENDDO
    142142
    143      print *, "Downloading data for tsurf_slope done"
     143     write(*,*) "Downloading data for tsurf_slope done"
    144144
    145145#ifndef CPP_STD
    … …  
    147147     if(soil_pem) then
    148148
    149      print *, "Downloading data for soiltemp_slope ..."
     149     write(*,*) "Downloading data for soiltemp_slope ..."
    150150
    151151DO islope=1,nslope
    … …  
    154154ENDDO
    155155
    156      print *, "Downloading data for soiltemp_slope done"
    157 
    158      print *, "Downloading data for watersoil_density ..."
     156     write(*,*) "Downloading data for soiltemp_slope done"
     157
     158     write(*,*) "Downloading data for watersoil_density ..."
    159159
    160160DO islope=1,nslope
    … …  
    163163ENDDO
    164164
    165      print *, "Downloading data for  watersoil_density  done"
    166 
    167      print *, "Downloading data for  watersurf_density  ..."
     165     write(*,*) "Downloading data for  watersoil_density  done"
     166
     167     write(*,*) "Downloading data for  watersurf_density  ..."
    168168
    169169DO islope=1,nslope
    … …  
    172172ENDDO
    173173
    174      print *, "Downloading data for  watersurf_density  done"
     174     write(*,*) "Downloading data for  watersurf_density  done"
    175175
    176176  endif !soil_pem
    … …  
    196196
    197197! Compute the minimum over the year for each point
    198   print *, "Computing the min of h2o_ice_slope"
     198  write(*,*) "Computing the min of h2o_ice_slope"
    199199  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4)
    200   print *, "Computing the min of co2_ice_slope"
     200  write(*,*) "Computing the min of co2_ice_slope"
    201201  min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4)
    202202
    203203!Compute averages
    204204
    205     print *, "Computing average of tsurf"
     205    write(*,*) "Computing average of tsurf"
    206206    tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
    207207
    … …  
    215215
    216216  if(soil_pem) then
    217     print *, "Computing average of tsoil"
     217    write(*,*) "Computing average of tsoil"
    218218    tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    219     print *, "Computing average of watersurf_density"
     219    write(*,*) "Computing average of watersurf_density"
    220220    watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen
    221221  endif

  • trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90

    r2980 r3050  
    101101do i=1,nvars
    102102  status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims,natts = numatts)
    103       print *, "namevar00= ", namevar
     103  write(*,*) "namevar00= ", namevar
    104104  if(status /= nf90_noerr) call handle_err(status)
    105105  allocate(dimid_var(numdims))
Note: See TracChangeset for help on using the changeset viewer.