Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3740) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3741) @@ -7,5 +7,5 @@ SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref, & nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & - play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,q_def_h2o_vap) + play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap) use ioipsl_getincom, only: getin ! To use 'getin' @@ -81,6 +81,7 @@ real, dimension(nlayer + 1), intent(out) :: plev ! intermediate pressure levels (pa) real, dimension(1), intent(out) :: latitude, longitude, cell_area -real, intent(out) :: atm_wat_profile, atm_wat_tau ! Force atmospheric water profiles -real, dimension(nlayer), intent(out) :: q_def_h2o_vap ! User-defined atmospheric water profile +logical, intent(out) :: prescribed_h2ovap ! Prescribe atmospheric water profile +real, intent(out) :: h2ovap_relax_time ! Relaxtion time towards the atmospheric water profile +real, dimension(nlayer), intent(out) :: q_prescribed_h2o_vap ! User-defined atmospheric water profile !======================================================================= @@ -778,59 +779,42 @@ if (.not. therestartfi .and. paleoclimate) perennial_co2ice = 0. -! Check if the atmospheric water profile is specified -! --------------------------------------------------- -! Adding an option to force atmospheric water values -atm_wat_profile = -1. ! Default: free atmospheric water profile -q_def_h2o_vap = 0. ! Default +! Prescription of atmospheric water profile +! ----------------------------------------- +prescribed_h2ovap = .false. ! Default: free atmospheric water profile +h2ovap_relax_time = -1. ! Default: no time relaxation +q_prescribed_h2o_vap = 0. ! Default if (water) then - write(*,*)'Force atmospheric water vapor profile (mixing ratio in kg/kg)?' - call getin('atm_wat_profile',atm_wat_profile) - write(*,*) 'atm_wat_profile = ', atm_wat_profile - if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1. + write(*,*)'Force atmospheric water profile (mixing ratio in kg/kg)?' + call getin('prescribed_h2ovap',prescribed_h2ovap) + write(*,*) 'prescribed_h2ovap = ', prescribed_h2ovap + if (prescribed_h2ovap) then + write(*,*) 'Atmospheric water profile is prescribed' + write(*,*) 'Unless it reaches saturation (maximal value)' + open(10,file = 'profile_prescribed_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr) + if (ierr == 0) then + do ilayer = 1,nlayer + read(10,*) q_prescribed_h2o_vap(ilayer) + enddo + else + error stop 'File "profile_prescribed_h2o_vap" was not found!' + endif + close(10) + + write(*,*) 'Relax atmospheric water profile with time constant (s)?' + + call getin('h2ovap_relax_time',h2ovap_relax_time) + write(*,*) 'h2ovap_relax_time = ', h2ovap_relax_time + if (h2ovap_relax_time < 0.) then + write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).' + else + write(*,*) 'Atmospheric water profile is relaxed towards the profile in "profile_prescribed_h2o_vap"' + write(*,*) 'Unless it reaches saturation (maximal value)' + endif + else write(*,*) 'Free atmospheric water vapor profile' write(*,*) 'Total water is conserved in the column' - else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0. - write(*,*) 'Dry atmospheric water vapor profile' - else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then ! if 0. <= and <= 1. - write(*,*) 'Prescribed uniform atmospheric water vapor profile' - write(*,*) 'Unless it reaches saturation (maximal value)' - else if (abs(atm_wat_profile + 2.) < 1.e-15) then ! if == -2. - write(*,*) 'Prescribed user-defined atmospheric water vapor profile' - write(*,*) 'Unless it reaches saturation (maximal value)' - open(10,file = 'profile_def_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr) - if (ierr == 0) then - do ilayer = 1,nlayer - read(10,*) q_def_h2o_vap(ilayer) - enddo - else - error stop 'File "profile_def_h2o_vap" was not found!' - endif - close(10) - else - error stop 'Water vapor profile value is not correct!' endif endif -! Check if the atmospheric water profile relaxation is specified -! -------------------------------------------------------------- -! Adding an option to relax atmospheric water values -atm_wat_tau = -1. ! Default: no time relaxation -if (water) then - write(*,*) 'Relax atmospheric water vapor profile with time constant (s)?' - call getin('atm_wat_tau',atm_wat_tau) - write(*,*) 'atm_wat_tau = ', atm_wat_tau - if (atm_wat_tau < 0.) then - write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).' - else - if (0. <= atm_wat_profile .and. atm_wat_profile <= 1. .or. abs(atm_wat_profile + 2.) < 1.e-15) then - write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile - write(*,*) 'Unless it reaches saturation (maximal value)' - else - write(*,*) 'Reference atmospheric water vapor profile not known!' - error stop 'Please, specify atm_wat_profile' - endif - endif -endif - END SUBROUTINE init_testphys1d Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3740) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3741) @@ -84,7 +84,8 @@ logical :: startfiles_1D, therestart1D, therestartfi, there -! Force atmospheric water profiles -real :: atm_wat_profile, atm_wat_tau -real, dimension(nlayer) :: q_def_h2o_vap +! Prescription of atmospheric water profile +logical :: prescribed_h2ovap +real :: h2ovap_relax_time +real, dimension(nlayer) :: q_prescribed_h2o_vap real, dimension(:), allocatable :: zqsat !======================================================================= @@ -150,5 +151,5 @@ call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, & time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & - play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,q_def_h2o_vap) + play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap) ! Write a "startfi" file @@ -196,20 +197,13 @@ enddo - ! Modify atmospheric water if asked - ! Added "atm_wat_profile" flag - ! ---------------------------- + ! Prescription of atmospheric water if asked + ! ------------------------------------------ if (water) then call watersat(nlayer,temp,play,zqsat) - if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then ! If atmospheric water is monitored with a uniform value - if (atm_wat_tau < 0.) then ! Case for prescribed atm_wat_profile: wet if >0, dry if =0 - q(1,:,igcm_h2o_vap) = min(zqsat,atm_wat_profile) - else ! Case for relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau - q(1,:,igcm_h2o_vap) = atm_wat_profile + (q(1,:,igcm_h2o_vap) - atm_wat_profile)*dexp(-dttestphys/atm_wat_tau) - endif - else if (abs(atm_wat_profile + 2.) < 1.e-15) then ! If atmospheric water is monitored with a user-defined water profile (with file "profile_def_h2o_vap") - if (atm_wat_tau < 0.) then ! Case for prescribed atm_wat_profile: wet if >0, dry if =0 - q(1,:,igcm_h2o_vap) = min(zqsat(:),q_def_h2o_vap(:)) - else ! Case for relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau - q(1,:,igcm_h2o_vap) = q_def_h2o_vap(:) + (q(1,:,igcm_h2o_vap) - q_def_h2o_vap(:))*dexp(-dttestphys/atm_wat_tau) + if (prescribed_h2ovap) then ! If atmospheric water profile is prescribed + if (h2ovap_relax_time < 0.) then ! Forcing + q(1,:,igcm_h2o_vap) = min(zqsat(:),q_prescribed_h2o_vap(:)) + else ! Relaxation + q(1,:,igcm_h2o_vap) = q_prescribed_h2o_vap(:) + (q(1,:,igcm_h2o_vap) - q_prescribed_h2o_vap(:))*dexp(-dttestphys/h2ovap_relax_time) endif endif