module read_profile_mod implicit none contains !=====================================================================================================================! ! SUBROUTINE: read_profile ===========================================================================================! !=====================================================================================================================! ! Author: Christophe Mathe ! Date: 09/06/2020 !---------------------------------------------------------------------------------------------------------------------! ! Subject: !--------- ! Read input profile of tracers listed in traceur.def !---------------------------------------------------------------------------------------------------------------------! ! Comments: !---------- ! - If no input profile found, q(traceur) and qsurf(traceur) set to 0, except: ! - q(co2) = 0.95 ! - q(hdo) = q(h2o)*2*155.76e-6*5 ! - q(dust_mass) and q(dust_number) initialized with Conrath dust ! ! - igcm is not available at this part of the code ! ! - To ensure that major isotopologue input profile is read before compute minor isotopologue profile, the ! calculation is performed at the end of the subroutine !---------------------------------------------------------------------------------------------------------------------! ! Algorithm: !----------- ! 1. Initialization of q and qsurf to 0 ! 2. Get indices of some tracers ! 3. Main ! 3.1. Try to open the input profile ! 3.1.1. Succeed ! 3.1.2. Fail ! 3.1.2.a. Some cases require a special initialization ! 4. Traitment for minor isotopologue if the major isotopologue input profile does not exist !=====================================================================================================================! subroutine read_profile(nb_tracer, nb_layer, qsurf, q, play) use infotrac, only: tname use tracer_mod, only: igcm_co2, igcm_h2o_vap, igcm_h2o_ice, & igcm_dust_number, igcm_dust_mass, & igcm_ccn_number, igcm_ccn_mass use aeropacity_mod, only: topdustref use dust_param_mod, only: odpref use comcstfi_h, only: pi implicit none !---------------------------------------------------------------------------------------------------------------------! ! VARIABLE DECLARATION !---------------------------------------------------------------------------------------------------------------------! ! Input arguments: !------------------ integer, intent(in) :: & nb_layer, & ! number of layer nb_tracer ! number of traceur read from traceur.def real, dimension(nb_layer), intent(in) :: play ! Pressure at the middle of the layers (Pa) !---------------------------------------------------------------------------------------------------------------------! ! Output arguments: !------------------- real, intent(out) :: & qsurf(nb_tracer), & ! kg/m2 q(nb_layer, nb_tracer) ! kg/kg of atmosphere !---------------------------------------------------------------------------------------------------------------------! ! Local: !-------- integer :: & iq, & ! loop over nb_tracer ilayer, & ! loop over nb_layer ierr, & ! open file iostat indice_h2o_vap = 0, & ! indice of h2o_vap tracer indice_h2o_ice = 0, & ! indice of h2o_ice tracer indice_hdo_vap = 0, & ! indice of hdo_vap tracer indice_hdo_ice = 0, & ! indice of hdo_ice tracer indice_dust_mass = 0, & ! indice of hdo_vap tracer indice_dust_number = 0 ! indice of hdo_ice tracer character(len=80), dimension(nb_tracer) :: & name_tracer ! array of all tracers already read in traceur.def logical :: & hdo_vap = .false., & ! used to compute hdo_vap profile if its input profile is missing (= .true.) hdo_ice = .false., & ! used to compute hdo_ice profile if its input profile is missing (= .true.) dust = .false. ! used to compute dust mass & number profiles if their input profiles are missing (= .true.) real :: zp, varian, r3n_q, ref_r0, r0_lift, mass2number_lift ! Dust variables real, parameter :: & nueff_lift = 0.5, & reff_lift = 3.e-6, & ! Effective radius of lifted dust (m) rho_dust = 2500 ! Mars dust density (kg.m-3) !=====================================================================================================================! !=== BEGIN ! !=====================================================================================================================! ! 1. Initialization of q and qsurf to 0 !---------------------------------------------------------------------------------------------------------------------! q(1:nb_layer,1:nb_tracer) = 0. qsurf(1:nb_tracer) = 0. name_tracer(1:nb_tracer) = "" !---------------------------------------------------------------------------------------------------------------------! ! 2. Get indices of some tracers: igcm is not available at this part of the code !---------------------------------------------------------------------------------------------------------------------! do iq = 1, nb_tracer write(name_tracer(iq),"(a)") tname(iq) if (trim(name_tracer(iq)) == 'h2o_vap') then indice_h2o_vap = iq else if (trim(name_tracer(iq)) == 'h2o_ice') then indice_h2o_ice = iq else if (trim(name_tracer(iq)) == 'hdo_vap') then indice_hdo_vap = iq else if (trim(name_tracer(iq)) == 'hdo_ice') then indice_hdo_ice = iq else if (trim(name_tracer(iq)) == 'dust_mass') then indice_dust_mass = iq else if (trim(name_tracer(iq)) == 'dust_number') then indice_dust_number = iq end if end do !---------------------------------------------------------------------------------------------------------------------! ! 3. Main !---------------------------------------------------------------------------------------------------------------------! do iq = 1, nb_tracer write(*,*)" tracer:",trim(name_tracer(iq)) !---------------------------------------------------------------------------------------------------------------------! ! 3.1. Try to open the input profile !---------------------------------------------------------------------------------------------------------------------! open(91,file='profile_'//trim(name_tracer(iq)), status='old', form='formatted', iostat=ierr) !---------------------------------------------------------------------------------------------------------------------! ! 3.1.1. Succeed: the input file exists !---------------------------------------------------------------------------------------------------------------------! if (ierr.eq.0) then read(91,*)qsurf(iq) if (trim(name_tracer(iq)) == 'co2') igcm_co2 = iq if (trim(name_tracer(iq)) == 'h2o_vap') igcm_h2o_vap = iq if (trim(name_tracer(iq)) == 'h2o_ice') igcm_h2o_ice = iq if (trim(name_tracer(iq)) == 'dust_number') igcm_dust_number = iq if (trim(name_tracer(iq)) == 'dust_mass') igcm_dust_mass = iq if (trim(name_tracer(iq)) == 'ccn_number') igcm_ccn_number = iq if (trim(name_tracer(iq)) == 'ccn_mass') igcm_ccn_mass = iq do ilayer = 1, nb_layer read(91,*)q(ilayer,iq) end do !---------------------------------------------------------------------------------------------------------------------! ! 3.1.2. Fail: the input file does not exist !---------------------------------------------------------------------------------------------------------------------! else write(*,*)"Warning: file profile_"//trim(name_tracer(iq))//" not found!" !---------------------------------------------------------------------------------------------------------------------! ! 3.1.2.a. Some cases require a special initialization !---------------------------------------------------------------------------------------------------------------------! select case(trim(name_tracer(iq))) case("co2") q(1:nb_layer,iq) = 0.95 write(*,*)"Initilization: q(co2) = 0.95" case("hdo_vap") hdo_vap = .true. write(*,*)"Initialization: q(hdo) = q(h2o)*2*155.76e-6*5" case("hdo_ice") hdo_ice = .true. write(*,*)"Initialization: q(hdo) = q(h2o)*2*155.76e-6*5" case("dust_mass","dust_number") dust = .true. write(*,*)"q(dust_mass) and q(dust_number) initialized with Conrath dust" case default write(*,*)"q("//trim(name_tracer(iq))//") and qsurf("//trim(name_tracer(iq))//") are set to 0!" end select end if close(91) end do ! of do iq = 1, nq !---------------------------------------------------------------------------------------------------------------------! ! 4. Traitment for minor isotopologue if the major isotopologue input profile does not exist. At this part, ensure that ! the main isotopologue profile is already read before minors isotopologues !---------------------------------------------------------------------------------------------------------------------! if (hdo_vap .and. indice_h2o_vap /= 0) then do ilayer = 1, nb_layer q(ilayer, indice_hdo_vap) = q(ilayer, indice_h2o_vap)*2*155.76e-6*5 end do end if if (hdo_ice .and. indice_h2o_ice /= 0) then qsurf(indice_hdo_ice) = qsurf(indice_h2o_ice) * 2*155.76e-6*5 do ilayer = 1, nb_layer q(ilayer, indice_hdo_ice) = q(ilayer, indice_h2o_ice) * 2*155.76e-6*5 end do end if if (dust) then varian = sqrt(log(1. + nueff_lift)) r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust) ref_r0 = exp(2.5*varian**2) r0_lift = reff_lift/ref_r0 mass2number_lift = r3n_q/r0_lift**3 do ilayer = 1,nb_layer zp = (odpref/play(ilayer))**(70./topdustref) q(ilayer,indice_dust_mass) = max(exp(0.007*(1. - max(zp,1.))),1.e-3) q(ilayer,indice_dust_number) = q(ilayer,indice_dust_mass)*mass2number_lift enddo endif !=====================================================================================================================! !=== END ! !=====================================================================================================================! end subroutine read_profile end module read_profile_mod