MODULE icolmdz_etat0 CONTAINS SUBROUTINE init_etat0 USE xios_mod USE omp_para USE getin_mod, ONLY: getin IMPLICIT NONE CHARACTER(LEN=256) :: etat0_lmdz ! @getin_doc etat0_database_type realm=initial type=String default=`legacy` ! When `etat0=="database"`, selects type of external file(s) from which initial state is read. Valid values : ! * `legacy` : reads from XIOS group `read_files_legacy`, see `xml/context_dynamico.xml` ! * `ERA5_forcing` : reads from XIOS group `read_files_ERA5_forcing`, see `xml/context_dynamico.xml` etat0_lmdz = "legacy" CALL getin("etat0_lmdz", etat0_lmdz) IF (is_omp_master) THEN CALL xios_set_fieldgroup_attr("read_fields", enabled=.TRUE.) IF (TRIM(etat0_lmdz) == "legacy") THEN CALL xios_set_filegroup_attr("read_files_legacy", enabled=.TRUE.) CALL xios_set_field_attr("relief_db", field_ref="relief_db_legacy") CALL xios_set_field_attr("ps_db", field_ref="ps_db_legacy") CALL xios_set_field_attr("z_db", field_ref="z_db_legacy") CALL xios_set_field_attr("ts_db", field_ref="ts_db_legacy") CALL xios_set_field_attr("u_db", field_ref="u_db_legacy") CALL xios_set_field_attr("v_db", field_ref="v_db_legacy") CALL xios_set_field_attr("temp_db", field_ref="temp_db_legacy") CALL xios_set_field_attr("q_db", field_ref="q_db_legacy") ELSE IF (TRIM(etat0_lmdz) == "ERA5_forcing") THEN CALL xios_set_filegroup_attr("read_files_ERA5_forcing", enabled=.TRUE.) CALL xios_set_field_attr("relief_db", field_ref="relief_db_forcing") CALL xios_set_field_attr("ps_db", field_ref="ps_db_forcing") CALL xios_set_field_attr("z_db", field_ref="z_db_forcing") CALL xios_set_field_attr("ts_db", field_ref="ts_db_forcing") CALL xios_set_field_attr("u_db", field_ref="u_db_forcing") CALL xios_set_field_attr("v_db", field_ref="v_db_forcing") CALL xios_set_field_attr("temp_db", field_ref="temp_db_forcing") CALL xios_set_field_attr("q_db", field_ref="q_db_forcing") ELSE PRINT*,"Bad selector for variable ",TRIM(etat0_lmdz), " option are (default), " STOP END IF END IF END SUBROUTINE init_etat0 SUBROUTINE etat0(f_ps, f_phis, f_theta_rhodz, f_u, f_q) USE icosa USE restart_mod USE wind_from_lonlat_mod USE write_field_mod USE time_mod USE transfert_mod USE xios_mod USE write_field_mod USE vertical_remap_mod USE theta2theta_rhodz_mod USE qsat_mod USE pression_mod USE omp_para USE tracer_icosa_mod IMPLICIT NONE TYPE(t_field), POINTER :: f_ps(:) TYPE(t_field), POINTER :: f_phis(:) TYPE(t_field), POINTER :: f_theta_rhodz(:) TYPE(t_field), POINTER :: f_u(:) TYPE(t_field), POINTER :: f_q(:) TYPE(t_field), POINTER, SAVE :: f_ulon_reg(:) TYPE(t_field), POINTER, SAVE :: f_ulat_reg(:) TYPE(t_field), POINTER, SAVE :: f_temp_reg(:) TYPE(t_field), POINTER, SAVE :: f_q_reg(:) TYPE(t_field), POINTER, SAVE :: f_ts(:) TYPE(t_field), POINTER, SAVE :: f_z(:) TYPE(t_field), POINTER, SAVE :: f_ulon(:) TYPE(t_field), POINTER, SAVE :: f_ulat(:) TYPE(t_field), POINTER, SAVE :: f_temp(:) TYPE(t_field), POINTER, SAVE :: f_q1(:) TYPE(t_field), POINTER, SAVE :: f_qsat(:) TYPE(t_field), POINTER, SAVE :: f_p(:) INTEGER :: nb_level REAL, ALLOCATABLE:: levels(:) INTEGER :: ind INTEGER :: iq, iq_H2Og, iq_H2Ol, iq_H2Os CALL xios_read_field("relief_db", f_phis) CALL writeField("relief_out", f_phis, once=.TRUE.) IF (is_omp_level_master) THEN DO ind = 1, ndomain IF (.NOT. assigned_domain(ind)) CYCLE f_phis(ind)%rval2d(:) = f_phis(ind)%rval2d(:)*g END DO END IF !$OMP BARRIER IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", n_glo=nb_level) CALL bcast_omp(nb_level) ALLOCATE (levels(nb_level)) IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", value=levels) CALL bcast_omp(levels) levels = levels*100 ! hectoPascal -> Pascal CALL allocate_field(f_ts, field_t, type_real, name="ts") CALL allocate_field(f_z, field_t, type_real, name="z") CALL allocate_field(f_ulon_reg, field_t, type_real, nb_level) CALL allocate_field(f_ulat_reg, field_t, type_real, nb_level) CALL allocate_field(f_temp_reg, field_t, type_real, nb_level) CALL allocate_field(f_q_reg, field_t, type_real, nb_level) CALL allocate_field(f_q1, field_t, type_real, llm) CALL allocate_field(f_qsat, field_t, type_real, llm) CALL allocate_field(f_p, field_t, type_real, llm + 1) CALL allocate_field(f_temp, field_t, type_real, llm) CALL allocate_field(f_ulon, field_t, type_real, llm) CALL allocate_field(f_ulat, field_t, type_real, llm) CALL xios_read_field("z_db", f_z) CALL xios_read_field("ps_db", f_ps) CALL xios_read_field("ts_db", f_ts) CALL writeField("ps_out", f_ps) !$OMP BARRIER ! CALL writeField("phis_out",f_phis,once=.TRUE.) ! CALL writeField("ts_out",f_ts,once=.TRUE.) ! make correction to ps due to relief at higher resolution ! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY IF (is_omp_level_master) THEN DO ind = 1, ndomain IF (.NOT. assigned_domain(ind)) CYCLE f_ps(ind)%rval2d(:) = f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:) - f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:)) END DO END IF !$OMP BARRIER CALL transfert_request(f_ps, req_i0) CALL writeField("ps_out", f_ps) CALL xios_read_field("temp_db", f_temp_reg) CALL vertical_remap(levels, f_temp_reg, f_ps, f_temp) CALL transfert_request(f_temp, req_i0) CALL temperature2theta_rhodz(f_ps, f_temp, f_theta_rhodz) CALL xios_read_field("u_db", f_ulon_reg) CALL vertical_remap(levels, f_ulon_reg, f_ps, f_ulon) CALL xios_read_field("v_db", f_ulat_reg) CALL vertical_remap(levels, f_ulat_reg, f_ps, f_ulat) CALL transfert_request(f_ulat, req_i0) CALL transfert_request(f_ulon, req_i0) CALL ulonlat2un(f_ulon, f_ulat, f_u) CALL xios_read_field("q_db", f_q_reg) CALL vertical_remap(levels, f_q_reg, f_ps, f_q1) CALL pression(f_ps, f_p) ! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2 CALL qsat(f_temp, f_p, f_qsat) CALL transfert_request(f_qsat, req_i0) DO iq = 1, nqtot IF (TRIM(tracers(iq)%name) == 'H2O_g') iq_H2Og=iq IF (TRIM(tracers(iq)%name) == 'H2O_l') iq_H2Ol=iq IF (TRIM(tracers(iq)%name) == 'H2O_s') iq_H2Os=iq ENDDO DO ind = 1, ndomain IF (.NOT. assigned_domain(ind)) CYCLE IF (tracers(iq_H2Ol)%has_default_init_value) THEN ! f_q(ind)%rval4d(:, :, iq_H2Ol) = tracers(iq_H2Ol)%default_init_value f_q(ind)%rval4d(:, :, iq_H2Ol) = tracers(iq_H2Ol)%default_init_value ELSE f_q(ind)%rval4d(:, :, iq_H2Ol) = 1e-6 END IF f_q(ind)%rval4d(:, :, iq_H2Og) = f_q1(ind)%rval3d(:, :)*f_qsat(ind)%rval3d(:, :)*0.01 WHERE(f_q(ind)%rval4d(:, :, iq_H2Og)<1e-20) f_q(ind)%rval4d(:, :, iq_H2Og) = 1e-20 f_q(ind)%rval4d(:, :, iq_H2Os) = 1e-10 END DO tracers(iq_H2Og)%already_initialized = .TRUE. tracers(iq_H2Ol)%already_initialized = .TRUE. tracers(iq_H2Os)%already_initialized = .TRUE. PRINT*,'CALLING ISOTOPE INITIALIZATION' CALL etat0_iso(f_q) ! CALL switch_iso(f_q) CALL writeField("tempdb_out", f_temp_reg) CALL writeField("temp_out", f_temp) CALL deallocate_field(f_ts) CALL deallocate_field(f_z) CALL deallocate_field(f_ulon_reg) CALL deallocate_field(f_ulat_reg) CALL deallocate_field(f_temp_reg) CALL deallocate_field(f_q_reg) CALL deallocate_field(f_q1) CALL deallocate_field(f_qsat) CALL deallocate_field(f_p) CALL deallocate_field(f_temp) CALL deallocate_field(f_ulon) CALL deallocate_field(f_ulat) END SUBROUTINE etat0 SUBROUTINE switch_iso(f_q) USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ USE isotopes_routines_mod, ONLY : fractcalk_liq, calcul_kcin ! ==> LMDZ USE infotrac_phy, ONLY: isoName, niso ! ==> LMDZ USE icosa ! ==> DYNAMICO USE tracer_icosa_mod USE field_mod USE omp_para, ONLY: is_master IMPLICIT NONE TYPE(t_field), POINTER :: f_q(:) INTEGER :: ind, iq, iq_H2Og, iq_H216Og, iq_H218Og, iq_H2Ol DO iq = 1, nqtot IF (TRIM(tracers(iq)%name) == 'H2O_l') iq_H2Ol=iq IF (TRIM(tracers(iq)%name) == 'H2O_g') iq_H2Og=iq IF (TRIM(tracers(iq)%name) == 'H216O_g') iq_H216Og=iq IF (TRIM(tracers(iq)%name) == 'H218O_g') iq_H218Og=iq ENDDO DO ind = 1, ndomain IF (.NOT. assigned_domain(ind)) CYCLE f_q(ind)%rval4d(:, :, iq_H2Ol) = f_q(ind)%rval4d(:, :, iq_H2Og) f_q(ind)%rval4d(:, :, iq_H218Og) = f_q(ind)%rval4d(:, :, iq_H2Og) f_q(ind)%rval4d(:, :, iq_H2Og) = 1. ENDDO END SUBROUTINE switch_iso SUBROUTINE etat0_iso(f_q) USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ USE isotopes_routines_mod, ONLY : fractcalk_liq, calcul_kcin ! ==> LMDZ USE infotrac_phy, ONLY: isoName, niso ! ==> LMDZ USE icosa ! ==> DYNAMICO USE tracer_icosa_mod USE field_mod USE omp_para, ONLY: is_master USE compute_transport_mod, ONLY : q_minimum IMPLICIT NONE TYPE(t_field), POINTER :: f_q(:) INTEGER :: iq, iqchild, iq_wiso CHARACTER(LEN=256) :: wiso_name ! idealised conditions for Rayleigh distillation ! R = R0_evap(wiso) * (q / q0_evap)**(alpha_ideal(wiso) - 1) REAL,PARAMETER :: q0_evap = 20e-3 ! initial condition of water vapour (kg/kg) : for Rayleigh distillation REAL,PARAMETER :: rh0_evap = 0.7 ! initial condition of relative humidity (1) : for R0_evap REAL :: ws0_evap = 2.0 ! initial condition of wind speed (m/s) : for R0_evap REAL :: kcin(niso) ! kinetic coefficient at evaporation = f(ws0_evap) : for R0_evap REAL,PARAMETER :: temp_ideal = 283. ! mean atmospheric temperature used for fractionation (K) : for Rayleigh distillation ! isotope-dependent REAL :: R0_evap ! mean value of initial isotopic ratio of evaporation over ocean : for Rayleigh distillation REAL :: alpha_ideal ! mean value of fractionation coefficient : for Rayleigh distillation REAL :: Rocean ! mean value of initial isotopic ratio of ocean : for R0_evap (m/s) REAL :: Rdefault_val ! default value of initial isotopic ratio REAL :: kcin_ocean ! kinetic coefficient at evaporation = f(ws0_evap) : for R0_evap INTEGER :: ind PRINT*, "ISOTOPE INITIALIZED ? ", using_iso IF (.NOT. using_iso) RETURN ! overwrite initialisation_iso => start from scratch initialisation_iso = 1 CALL calcul_kcin(ws0_evap, kcin) DO ind = 1, ndomain IF (.NOT. assigned_domain(ind)) CYCLE ! CAa : loop on the tracers: look for water components DO iq = 1, nqtot ! looking for vapour = 'H2O_g' ! "from readTracFiles_mod: (Only convertable names are water descendants names H2O_, _, __, with:)" IF (TRIM(tracers(iq)%name) == 'H2O_g') THEN ! Loop on children of water vapour ! ================================ ! idealised Rayleigh distillation DO iqchild = 1, tracers(iq)%nb_children iq_wiso = tracers(iq)%children(iqchild) wiso_name = tracers(iq_wiso)%name SELECT CASE(TRIM(wiso_name)) CASE('H216O_g'); CALL fractcalk_liq(iso_eau, temp_ideal, alpha_ideal) Rocean = toce(iso_eau) kcin_ocean= kcin(iso_eau) CASE('H217O_g') CALL fractcalk_liq(iso_O17, temp_ideal, alpha_ideal) Rocean = toce(iso_O17) kcin_ocean= kcin(iso_O17) CASE('H218O_g') CALL fractcalk_liq(iso_O18, temp_ideal, alpha_ideal) Rocean = toce(iso_O18) kcin_ocean= kcin(iso_O18) CASE('HDO_g') CALL fractcalk_liq(iso_HDO, temp_ideal, alpha_ideal) Rocean = toce(iso_HDO) kcin_ocean= kcin(iso_O18) CASE('HTO_g') CALL fractcalk_liq(iso_HTO, temp_ideal, alpha_ideal) Rocean = toce(iso_HTO) kcin_ocean= kcin(iso_O18) CASE DEFAULT PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP' STOP END SELECT ! R0_evap = R evaporation from Merlivat R0_evap = Rocean / alpha_ideal * (1.0 - kcin_ocean) / (1.0 - kcin_ocean * rh0_evap) ! R_iso = R from idealised Rayleigh distillation ! R_iso = R0_evap(wiso) * (q / q0_evap)**(alpha_ideal(wiso) - 1) f_q(ind)%rval4d(:, :, iq_wiso) = R0_evap * (min(q0_evap, f_q(ind)%rval4d(:, :, iq)) / q0_evap)**(alpha_ideal - 1) IF (is_master) PRINT*,TRIM(wiso_name)," ",Rocean, kcin_ocean, R0_evap, alpha_ideal, MAXVAL(f_q(ind)%rval4d(:, :, iq)),MINVAL(f_q(ind)%rval4d(:, :, iq)), MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso)) tracers(iq_wiso)%already_initialized = .TRUE. END DO ELSE IF (TRIM(tracers(iq)%name) == 'H2O_l') THEN ! TODO : update, it's wrong DO iqchild = 1, tracers(iq)%nb_children iq_wiso = tracers(iq)%children(iqchild) wiso_name = tracers(iq_wiso)%name SELECT CASE(wiso_name) CASE('H216O_l') Rdefault_val = Rdefault(iso_eau) CASE('H217O_l') Rdefault_val = Rdefault(iso_O17) CASE('H218O_l') Rdefault_val = Rdefault(iso_O18) CASE('HDO_l') Rdefault_val = Rdefault(iso_HDO) CASE('HTO_l') Rdefault_val = Rdefault(iso_HTO) CASE DEFAULT PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP' STOP END SELECT f_q(ind)%rval4d(:, :, iq_wiso) = Rdefault_val IF (is_master) PRINT*,TRIM(wiso_name)," ",Rdefault_val, MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso)) tracers(iq_wiso)%already_initialized = .TRUE. END DO ELSE IF (TRIM(tracers(iq)%name) == 'H2O_s') THEN DO iqchild = 1, tracers(iq)%nb_children iq_wiso = tracers(iq)%children(iqchild) wiso_name = tracers(iq_wiso)%name SELECT CASE(wiso_name) CASE('H216O_s') Rdefault_val = Rdefault(iso_eau) CASE('H217O_s') Rdefault_val = Rdefault(iso_O17) CASE('H218O_s') Rdefault_val = Rdefault(iso_O18) CASE('HDO_s') Rdefault_val = Rdefault(iso_HDO) CASE('HTO_s') Rdefault_val = Rdefault(iso_HTO) CASE DEFAULT PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP' STOP END SELECT f_q(ind)%rval4d(:, :, iq_wiso) = Rdefault_val IF (is_master) PRINT*,TRIM(wiso_name)," ",Rdefault_val, MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso)) tracers(iq_wiso)%already_initialized = .TRUE. END DO END IF ENDDO ENDDO END SUBROUTINE etat0_iso END MODULE icolmdz_etat0