Index: ICOSA_LMDZ/src/phylmdiso/icolmdz_etat0.f90 =================================================================== --- ICOSA_LMDZ/src/phylmdiso/icolmdz_etat0.f90 (revision 5590) +++ ICOSA_LMDZ/src/phylmdiso/icolmdz_etat0.f90 (revision 5590) @@ -0,0 +1,396 @@ +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 Index: ICOSA_LMDZ/src/phylmdiso/icolmdz_param_gravity_wave.f90 =================================================================== --- ICOSA_LMDZ/src/phylmdiso/icolmdz_param_gravity_wave.f90 (revision 5590) +++ ICOSA_LMDZ/src/phylmdiso/icolmdz_param_gravity_wave.f90 (revision 5590) @@ -0,0 +1,593 @@ +MODULE icolmdz_param_gravity_wave + + INTEGER, PARAMETER :: gw_legacy=1 + INTEGER, PARAMETER :: gw_sso=2 + INTEGER,SAVE :: gw_method + !$OMP THREADPRIVATE(gw_method) + +CONTAINS + + SUBROUTINE init_param_gravity_wave + USE getin_mod + USE xios_mod + IMPLICIT NONE + CHARACTER(LEN=255) :: param_gw_method + LOGICAL :: create_etat0_limit + + create_etat0_limit = .FALSE. + CALL getin('create_etat0_limit', create_etat0_limit) + IF (.NOT. create_etat0_limit) RETURN + + param_gw_method='sso' + CALL getin('param_gw_method', param_gw_method) + SELECT CASE (TRIM(param_gw_method)) + CASE ('legacy') + gw_method = gw_legacy + CASE ('sso') + gw_method = gw_sso + CASE default + PRINT *, 'Bad selector for param_gw_method : <', TRIM(param_gw_method), & + ' > options are , ' + STOP + END SELECT + + SELECT CASE (gw_method) + CASE (gw_legacy) + !$OMP MASTER + CALL xios_set_file_attr("relief_gw",enabled=.TRUE.) + CALL xios_set_fieldgroup_attr("gw_read_access",read_access=.TRUE.) + !$OMP END MASTER + CASE (gw_sso) + !$OMP MASTER + CALL xios_set_file_attr("orography",enabled=.TRUE.) + CALL xios_set_fieldgroup_attr("gwsso_read_access",read_access=.TRUE.) + !$OMP END MASTER + END SELECT + + END SUBROUTINE init_param_gravity_wave + + SUBROUTINE param_gravity_wave + USE getin_mod + IMPLICIT NONE + LOGICAL :: create_etat0_limit + + create_etat0_limit = .FALSE. + CALL getin('create_etat0_limit', create_etat0_limit) + IF (.NOT. create_etat0_limit) RETURN + + SELECT CASE (gw_method) + CASE (gw_legacy) + CALL param_gravity_wave_legacy + CASE (gw_sso) + CALL param_gravity_wave_sso + END SELECT + + END SUBROUTINE param_gravity_wave + + SUBROUTINE param_gravity_wave_sso + ! from dynamico + USE icosa + USE xios_mod + USE mpipara + USE earth_const + USE transfert_mod + ! from lmdz physic + USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw + USE dimphy, ONLY : klon + !from icosa_lmdz + USE distrib_icosa_lmdz_mod + IMPLICIT NONE + TYPE(t_field), POINTER, SAVE :: f_zmea(:), f_zpic(:), f_zval(:), f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:) + REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:) + + CALL allocate_field(f_zmea, field_t, type_real) + CALL allocate_field(f_zpic, field_t, type_real) + CALL allocate_field(f_zval, field_t, type_real) + CALL allocate_field(f_zstd, field_t, type_real) + CALL allocate_field(f_zsig, field_t, type_real) + CALL allocate_field(f_zgam, field_t, type_real) + CALL allocate_field(f_zthe, field_t, type_real) + + CALL xios_read_field("zmea_sso",f_zmea) + CALL transfert_request(f_zmea, req_i0) + CALL transfert_request(f_zmea, req_i1) + + CALL xios_read_field("zstd_sso",f_zstd) + CALL transfert_request(f_zstd, req_i0) + CALL transfert_request(f_zstd, req_i1) + + CALL xios_read_field("zsig_sso",f_zsig) + CALL transfert_request(f_zsig, req_i0) + CALL transfert_request(f_zsig, req_i1) + + CALL xios_read_field("zgam_sso",f_zgam) + CALL transfert_request(f_zgam, req_i0) + CALL transfert_request(f_zgam, req_i1) + + CALL xios_read_field("zthe_sso",f_zthe) + CALL transfert_request(f_zthe, req_i0) + CALL transfert_request(f_zthe, req_i1) + + ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon)) + CALL transfer_icosa_to_lmdz(f_zmea , zmea_phy) + CALL transfer_icosa_to_lmdz(f_zstd , zstd_phy) + CALL transfer_icosa_to_lmdz(f_zsig , zsig_phy) + CALL transfer_icosa_to_lmdz(f_zgam , zgam_phy) + CALL transfer_icosa_to_lmdz(f_zthe , zthe_phy) + zpic_phy(:) = zmea_phy(:) + zstd_phy(:) + zval_phy(:) = MAX(zmea_phy(:) - zstd_phy(:), 0.) + CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy) + + END SUBROUTINE param_gravity_wave_sso + + SUBROUTINE param_gravity_wave_legacy + ! from dynamico + USE icosa + USE xios_mod + USE mpipara + USE earth_const + USE transfert_mod + ! from lmdz physic + USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw + USE dimphy, ONLY : klon + !from icosa_lmdz + USE distrib_icosa_lmdz_mod + IMPLICIT NONE + INTEGER :: ibegin, jbegin, ni, nj, ni_glo,nj_glo + REAL, ALLOCATABLE :: z(:,:) + REAL, ALLOCATABLE :: lon1d(:), lon2d_in(:,:), lon2d_out(:,:) + REAL, ALLOCATABLE :: lat1d(:), lat2d_in(:,:), lat2d_out(:,:) + REAL,ALLOCATABLE :: delta_x(:), delta_y(:) + REAL, ALLOCATABLE :: mask_in(:,:) + REAL, ALLOCATABLE :: zmea_in(:,:), zpic_in(:,:), zval_in(:,:), ztz_in(:,:), zytzy_in(:,:), zxtzx_in(:,:), zxtzy_in(:,:) + + TYPE(t_field), POINTER, SAVE :: f_mask(:), f_zmea(:), f_zpic(:), f_zval(:), f_ztz(:), f_zytzy(:), f_zxtzx(:), f_zxtzy(:) + TYPE(t_field), POINTER, SAVE :: f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:) + REAL(rstd),POINTER :: mask(:), zmea(:), zpic(:), zval(:), ztz(:), zytzy(:), zxtzx(:), zxtzy(:), zstd(:), zsig(:), zgam(:), zthe(:) + REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:) + + INTEGER :: i, j, ind + + CALL allocate_field(f_mask, field_t, type_real) + CALL allocate_field(f_zmea, field_t, type_real) + CALL allocate_field(f_zpic, field_t, type_real) + CALL allocate_field(f_zval, field_t, type_real) + CALL allocate_field(f_ztz, field_t, type_real) + CALL allocate_field(f_zytzy, field_t, type_real) + CALL allocate_field(f_zxtzx, field_t, type_real) + CALL allocate_field(f_zxtzy, field_t, type_real) + CALL allocate_field(f_zstd, field_t, type_real) + CALL allocate_field(f_zsig, field_t, type_real) + CALL allocate_field(f_zgam, field_t, type_real) + CALL allocate_field(f_zthe, field_t, type_real) + +!$OMP MASTER + CALL xios_get_domain_attr("domain_relief_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj, nj_glo=nj_glo) + PRINT*, ni,nj,ni_glo,nj_glo + + ALLOCATE( lon1d(ni), lat1d(nj)) + CALL xios_get_domain_attr("domain_relief_gw",lonvalue_1d=lon1d, latvalue_1d=lat1d) + + ALLOCATE(z(0:ni+1,0:nj+1)) + CALL xios_recv_field("relief_exp", z) + + ALLOCATE(lon2d_in(ni,nj), lat2d_in(ni,nj) ) + ALLOCATE(lon2d_out(0:ni+1,0:nj+1), lat2d_out(0:ni+1,0:nj+1) ) + + DO i=1,ni + lon2d_in(i,:) = lon1d(i) + END DO + + DO j=1,nj + lat2d_in(:,j) = lat1d(j) + END DO + + CALL xios_send_field("lon2d_in", lon2d_in) + CALL xios_recv_field("lon2d_out", lon2d_out) + CALL xios_send_field("lat2d_in", lat2d_in) + CALL xios_recv_field("lat2d_out", lat2d_out) + + ALLOCATE(delta_x(nj)) + ALLOCATE(delta_y(nj)) + + delta_y(:) = Pi*radius / nj_glo !! difference with lmdz reference : delta_y=2*Pi*radius/nj_glo ?? + + DO j=1,nj + delta_x(j) = 2*Pi*radius / ni_glo * cos( lat2d_out(1,j) * pi / 180.) + ENDDO + + !north pole + IF (jbegin==0) THEN + z(:,0) = z(:,1) + delta_y(1) = delta_y(1)*0.5 + ENDIF + + !south pole + IF (jbegin+nj==nj_glo) THEN + z(:,nj+1) = z(:,nj) + delta_y(nj) = delta_y(nj)*0.5 + ENDIF + + ALLOCATE(mask_in(ni,nj)) + ALLOCATE(zmea_in(ni,nj), zpic_in(ni,nj), zval_in(ni,nj), ztz_in(ni,nj), zytzy_in(ni,nj), zxtzx_in(ni,nj), zxtzy_in(ni,nj) ) + + mask_in(:,:)=0 + DO j=1,nj + DO i=1,ni + IF ( z(i,j)>1. ) mask_in(i,j) = 1. + zmea_in(i,j) = z(i,j) + zpic_in(i,j) = z(i,j) + zval_in(i,j) = z(i,j) + ztz_in(i,j) = z(i,j)*z(i,j) + zytzy_in(i,j) = (z(i,j+1)-z(i,j-1))**2/(2*delta_y(j))**2 + zxtzx_in(i,j) = (z(i+1,j)-z(i-1,j))**2/(2*delta_x(j))**2 + zxtzy_in(i,j) = (z(i,j+1)-z(i,j-1)) / (2*delta_y(j)) *(z(i+1,j)-z(i-1,j)) / (2*delta_x(j)) + ENDDO + ENDDO + + CALL xios_send_field("mask_in",mask_in) + CALL xios_send_field("zmea_in",zmea_in) + CALL xios_send_field("zpic_in",zpic_in) + CALL xios_send_field("zval_in",zval_in) + CALL xios_send_field("ztz_in",ztz_in) + CALL xios_send_field("zytzy_in",zytzy_in) + CALL xios_send_field("zxtzx_in",zxtzx_in) + CALL xios_send_field("zxtzy_in",zxtzy_in) + + !! only for test, to compare projection on ico grid Vs lmdz_reg grid + !! CALL compute_regular_param + + !$OMP END MASTER + + CALL xios_read_field("mask_out", f_mask) + CALL transfert_request(f_mask, req_i0) + CALL transfert_request(f_mask, req_i1) + + CALL xios_read_field("zmea_out",f_zmea) + CALL transfert_request(f_zmea, req_i0) + CALL transfert_request(f_zmea, req_i1) + + CALL xios_read_field("zpic_out",f_zpic) + CALL transfert_request(f_zpic, req_i0) + CALL transfert_request(f_zpic, req_i1) + + CALL xios_read_field("zval_out",f_zval) + CALL transfert_request(f_zval, req_i0) + CALL transfert_request(f_zval, req_i1) + + CALL xios_read_field("ztz_out", f_ztz) + CALL transfert_request(f_ztz, req_i0) + CALL transfert_request(f_ztz, req_i1) + + CALL xios_read_field("zytzy_out", f_zytzy) + CALL transfert_request(f_zytzy, req_i0) + CALL transfert_request(f_zytzy, req_i1) + + CALL xios_read_field("zxtzx_out", f_zxtzx) + CALL transfert_request(f_zxtzx, req_i0) + CALL transfert_request(f_zxtzx, req_i1) + + CALL xios_read_field("zxtzy_out", f_zxtzy) + CALL transfert_request(f_zxtzy, req_i0) + CALL transfert_request(f_zxtzy, req_i1) + + !$OMP BARRIER + !$OMP MASTER + + !! compute zstd + + DO ind=1,ndomain + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + zmea =f_zmea(ind) + zstd = f_zstd(ind) + ztz = f_ztz(ind) + CALL compute_zstd(zmea, ztz, zstd) + ENDDO + + !! smmothing ? + + DO ind=1,ndomain + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + + zmea=f_zmea(ind) + zpic=f_zpic(ind) + zval=f_zval(ind) + zstd=f_zstd(ind) + zytzy=f_zytzy(ind) + zxtzx=f_zxtzx(ind) + zxtzy=f_zxtzy(ind) + + CALL compute_smoothing(zmea) + CALL compute_smoothing(zpic) + CALL compute_smoothing(zval) + CALL compute_smoothing(zstd) + CALL compute_smoothing(zytzy) + CALL compute_smoothing(zxtzx) + CALL compute_smoothing(zxtzy) + + ENDDO + + !! SLOPE, ANISOTROPY AND THETA ANGLE + + DO ind=1,ndomain + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + + zytzy=f_zytzy(ind) + zxtzx=f_zxtzx(ind) + zxtzy=f_zxtzy(ind) + zsig=f_zsig(ind) + zgam=f_zgam(ind) + zthe=f_zthe(ind) + CALL compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe) + + ENDDO + + !! apply mask + + DO ind=1,ndomain + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + + mask=f_mask(ind) + zmea=f_zmea(ind) + zpic=f_zpic(ind) + zval=f_zval(ind) + zstd=f_zstd(ind) + zsig=f_zsig(ind) + zgam=f_zgam(ind) + zthe=f_zthe(ind) + + CALL compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe) + + ENDDO + +!$OMP END MASTER +!$OMP BARRIER + + !! for debugging + !! CALL xios_write_field("zmask", f_mask) + !! CALL xios_write_field("zmea", f_zmea) + !! CALL xios_write_field("zpic", f_zpic) + !! CALL xios_write_field("zval", f_zval) + !! CALL xios_write_field("zsig", f_zsig) + !! CALL xios_write_field("zgam", f_zgam) + !! CALL xios_write_field("zthe", f_zthe) + + ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon)) + CALL transfer_icosa_to_lmdz(f_zmea , zmea_phy) + CALL transfer_icosa_to_lmdz(f_zpic , zpic_phy) + CALL transfer_icosa_to_lmdz(f_zval , zval_phy) + CALL transfer_icosa_to_lmdz(f_zstd , zstd_phy) + CALL transfer_icosa_to_lmdz(f_zsig , zsig_phy) + CALL transfer_icosa_to_lmdz(f_zgam , zgam_phy) + CALL transfer_icosa_to_lmdz(f_zthe , zthe_phy) + CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy) + + END SUBROUTINE param_gravity_wave_legacy + + SUBROUTINE compute_smoothing(X) + USE icosa + IMPLICIT NONE + REAL(rstd),INTENT(INOUT) :: X(iim*jjm) + REAL(rstd) :: tmp(iim*jjm) + INTEGER :: i,j,ij + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + tmp(ij) = (X(ij) + (X(ij+t_right) + X(ij+t_rup) + X(ij+t_lup) + X(ij+t_left) + X(ij+t_ldown) + X(ij+t_rdown))/3.)/3. + ENDDO + ENDDO + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + X(ij) = tmp(ij) + ENDDO + ENDDO + + END SUBROUTINE + + SUBROUTINE compute_zstd(zmea, ztz, zstd) + USE icosa + IMPLICIT NONE + REAL(rstd),INTENT(IN) :: zmea(iim*jjm) + REAL(rstd),INTENT(IN) :: ztz(iim*jjm) + REAL(rstd),INTENT(OUT) :: zstd(iim*jjm) + + INTEGER :: i,j,ij + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + zstd(ij) = ztz(ij) - zmea(ij)*zmea(ij) + IF (zstd(ij)<0) zstd(ij)=0 + zstd(ij)=SQRT(zstd(ij)) + ENDDO + ENDDO + + END SUBROUTINE compute_zstd + + + SUBROUTINE compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe) + USE icosa + IMPLICIT NONE + REAL(rstd),INTENT(IN) :: zxtzx(iim*jjm) + REAL(rstd),INTENT(IN) :: zytzy(iim*jjm) + REAL(rstd),INTENT(IN) :: zxtzy(iim*jjm) + REAL(rstd),INTENT(OUT) :: zsig(iim*jjm) + REAL(rstd),INTENT(OUT) :: zgam(iim*jjm) + REAL(rstd),INTENT(OUT) :: zthe(iim*jjm) + + REAL(rstd) :: xk,xl,xm,xp,xq,xw + + INTEGER :: i,j,ij + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + + ij=(j-1)*iim+i + xk=(zxtzx(ij)+zytzy(ij))/2. + xl=(zxtzx(ij)-zytzy(ij))/2. + xm=zxtzy(ij) + xp=xk-SQRT(xl**2+xm**2) + xq=xk+SQRT(xl**2+xm**2) + xw=1.e-8 + IF(xp<=xw) xp=0. + IF(xq<=xw) xq=xw + IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm) + + !--- SLOPE, ANISOTROPY AND THETA ANGLE + + zsig(ij)=SQRT(xq) + zgam(ij)=xp/xq + zthe(ij)=90.*ATAN2(xm,xl)/Pi + + ENDDO + ENDDO + + END SUBROUTINE compute_sigma_gamma_theta + + SUBROUTINE compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe) + USE icosa + IMPLICIT NONE + REAL(rstd),INTENT(IN) :: mask(iim*jjm) + REAL(rstd),INTENT(OUT) :: zmea(iim*jjm) + REAL(rstd),INTENT(OUT) :: zpic(iim*jjm) + REAL(rstd),INTENT(OUT) :: zval(iim*jjm) + REAL(rstd),INTENT(OUT) :: zstd(iim*jjm) + REAL(rstd),INTENT(OUT) :: zsig(iim*jjm) + REAL(rstd),INTENT(OUT) :: zgam(iim*jjm) + REAL(rstd),INTENT(OUT) :: zthe(iim*jjm) + + + INTEGER :: i,j,ij + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + IF (mask(ij)<0.1) THEN + zmea(ij)=0. + zpic(ij)=0. + zval(ij)=0. + zstd(ij)=0. + zsig(ij)=0. + zgam(ij)=0. + zthe(ij)=0. + ENDIF + ENDDO + ENDDO + + END SUBROUTINE compute_masking + + SUBROUTINE compute_regular_param + USE xios_mod + USE icosa + IMPLICIT NONE + REAL, ALLOCATABLE :: mask(:,:) , zmea(:,:), zpic(:,:), zval(:,:), ztz(:,:), zstd(:,:) + REAL, ALLOCATABLE :: zytzy(:,:) , zxtzx(:,:), zxtzy(:,:), zsig(:,:) , zgam(:,:), zthe(:,:) + REAL :: xk, xl, xm, xp, xq, xw + INTEGER :: ni, nj, ibegin, jbegin, ni_glo, nj_glo + INTEGER :: i,j + + CALL xios_get_domain_attr("regular_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj, nj_glo=nj_glo) + + ALLOCATE(mask(0:ni+1,0:nj+1), zmea(0:ni+1,0:nj+1), zpic(0:ni+1,0:nj+1), zval(0:ni+1,0:nj+1), ztz(0:ni+1,0:nj+1)) + ALLOCATE(zstd(0:ni+1,0:nj+1), zytzy(0:ni+1,0:nj+1), zxtzx(0:ni+1,0:nj+1), zxtzy(0:ni+1,0:nj+1)) + ALLOCATE(zsig(0:ni+1,0:nj+1), zgam(0:ni+1,0:nj+1), zthe(0:ni+1,0:nj+1)) + + CALL xios_recv_field("mask_reg_exp",mask) + CALL xios_recv_field("zmea_reg_exp",zmea) + CALL xios_recv_field("zpic_reg_exp",zpic) + CALL xios_recv_field("zval_reg_exp",zval) + CALL xios_recv_field("ztz_reg_exp",ztz) + CALL xios_recv_field("zytzy_reg_exp",zytzy) + CALL xios_recv_field("zxtzx_reg_exp",zxtzx) + CALL xios_recv_field("zxtzy_reg_exp",zxtzy) + + DO j=0,nj+1 + DO i=0, ni+1 + zstd(i,j) = ztz(i,j)-zmea(i,j)*zmea(i,j) + ENDDO + ENDDO + + ! smoothing + CALL smoothing_regular(ni,nj, zmea) + CALL smoothing_regular(ni,nj, zpic) + CALL smoothing_regular(ni,nj, zval) + CALL smoothing_regular(ni,nj, zstd) + CALL smoothing_regular(ni,nj, zytzy) + CALL smoothing_regular(ni,nj, zxtzx) + CALL smoothing_regular(ni,nj, zxtzy) + + DO j=1,nj + DO i=1, ni + xk=(zxtzx(i,j)+zytzy(i,j))/2. + xl=(zxtzx(i,j)-zytzy(i,j))/2. + xm=zxtzy(i,j) + xp=xk-SQRT(xl**2+xm**2) + xq=xk+SQRT(xl**2+xm**2) + xw=1.e-8 + IF(xp<=xw) xp=0. + IF(xq<=xw) xq=xw + IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm) + + !--- SLOPE, ANISOTROPY AND THETA ANGLE + + zsig(i,j)=SQRT(xq) + zgam(i,j)=xp/xq + zthe(i,j)=90.*ATAN2(xm,xl)/Pi + ENDDO + ENDDO + + DO j=1,nj + DO i=1,ni + IF (mask(i,j)<0.1) THEN + zmea(i,j)=0. + zpic(i,j)=0. + zval(i,j)=0. + zstd(i,j)=0. + zsig(i,j)=0. + zgam(i,j)=0. + zthe(i,j)=0. + ENDIF + ENDDO + ENDDO + + CALL xios_send_field("zmask_regout", mask(1:ni,1:nj)) + CALL xios_send_field("zmea_regout", zmea(1:ni,1:nj)) + CALL xios_send_field("zpic_regout", zpic(1:ni,1:nj)) + CALL xios_send_field("zval_regout", zval(1:ni,1:nj)) + CALL xios_send_field("zsig_regout", zsig(1:ni,1:nj)) + CALL xios_send_field("zgam_regout", zgam(1:ni,1:nj)) + CALL xios_send_field("zthe_regout", zthe(1:ni,1:nj)) + + END SUBROUTINE compute_regular_param + + SUBROUTINE smoothing_regular(ni,nj,var) + IMPLICIT NONE + INTEGER, INTENT(IN) :: ni, nj + REAL, INTENT(INOUT) :: var(0:ni+1,0:nj+1) + REAL :: tmp(0:ni+1,0:nj+1) + INTEGER :: i,j + + DO j=1,nj + DO i=1,ni + tmp(i,j) = (var(i,j) + 0.5*(var(i+1,j) + var(i-1,j) + var(i,j+1) +var(i,j-1)) & + + 0.25*(var(i+1,j+1) + var(i-1,j+1) + var(i+1,j-1) + var(i-1,j-1)))/4. + ENDDO + ENDDO + + DO j=1,nj + DO i=1,ni + var(i,j)=tmp(i,j) + ENDDO + ENDDO + + END SUBROUTINE smoothing_regular + + +END MODULE icolmdz_param_gravity_wave Index: ICOSA_LMDZ/src/phylmdiso/interface_icosa_lmdz.F90 =================================================================== --- ICOSA_LMDZ/src/phylmdiso/interface_icosa_lmdz.F90 (revision 5590) +++ ICOSA_LMDZ/src/phylmdiso/interface_icosa_lmdz.F90 (revision 5590) @@ -0,0 +1,1342 @@ +MODULE interface_icosa_lmdz_mod + + USE field_mod, ONLY: t_field + USE transfert_mod, ONLY: t_message + + + TYPE(t_message),SAVE :: req_u, req_z + TYPE(t_message),SAVE :: req_dps0, req_dulon0, req_dulat0, req_dTemp0, req_dq0 + + TYPE(t_field),POINTER,SAVE :: f_p(:) + TYPE(t_field),POINTER,SAVE :: f_pks(:) + TYPE(t_field),POINTER,SAVE :: f_pk(:) + TYPE(t_field),POINTER,SAVE :: f_p_layer(:) + TYPE(t_field),POINTER,SAVE :: f_theta(:) + TYPE(t_field),POINTER,SAVE :: f_phi(:) + TYPE(t_field),POINTER,SAVE :: f_Temp(:) + TYPE(t_field),POINTER,SAVE :: f_ulon(:) + TYPE(t_field),POINTER,SAVE :: f_ulat(:) + TYPE(t_field),POINTER,SAVE :: f_vort(:) + TYPE(t_field),POINTER,SAVE :: f_vortc(:) + TYPE(t_field),POINTER,SAVE :: f_dulon(:) + TYPE(t_field),POINTER,SAVE :: f_dulat(:) + TYPE(t_field),POINTER,SAVE :: f_dTemp(:) + TYPE(t_field),POINTER,SAVE :: f_dq(:) + TYPE(t_field),POINTER,SAVE :: f_dps(:) + TYPE(t_field),POINTER,SAVE :: f_duc(:) + TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:) + TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:) + + INTEGER :: start_clock + INTEGER :: stop_clock + INTEGER :: count_clock=0 + + INTEGER,SAVE :: nbp_phys + INTEGER,SAVE :: nbp_phys_glo + + REAL,ALLOCATABLE,SAVE :: q_ave(:, :) + !$OMP THREADPRIVATE(q_ave) + + REAL, PARAMETER :: q_epsilon = 1e-8 + +CONTAINS + + SUBROUTINE pre_initialize_physics + USE etat0_plugin_mod + USE icolmdz_etat0, ONLY: init_etat0_lmdz => init_etat0, etat0_lmdz => etat0 + USE icolmdz_param_gravity_wave, ONLY: init_param_gravity_wave + USE isotopes_mod, ONLY : using_iso + IMPLICIT NONE + !$OMP PARALLEL + + init_etat0_plugin => init_etat0_lmdz + etat0_plugin => etat0_lmdz + + CALL init_param_gravity_wave + !$OMP END PARALLEL + END SUBROUTINE pre_initialize_physics + + SUBROUTINE initialize_physics + USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz +! from dynamico + USE domain_mod + USE dimensions + USE mpi_mod + USE mpipara + USE disvert_mod + USE xios_mod + USE time_mod , init_time_icosa=> init_time + USE transfert_mod + USE nudging_mod, ONLY : lam_halo_scheme + +! from LMDZ + USE mod_grid_phy_lmdz, ONLY : unstructured + USE mod_phys_lmdz_omp_data, ONLY: klon_omp + USE transfert_mod + USE physics_distribution_mod, ONLY : init_physics_distribution + USE infotrac_phy, ONLY: init_infotrac_phy + USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave + + + + IMPLICIT NONE + INTEGER :: ind,i,j,ij,pos,h + REAL(rstd),POINTER :: bounds_lon(:,:) + REAL(rstd),POINTER :: bounds_lat(:,:) + + REAL(rstd),ALLOCATABLE :: latfi(:) + REAL(rstd),ALLOCATABLE :: lonfi(:) + REAL(rstd),ALLOCATABLE :: airefi(:) + REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:) + REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:) + LOGICAL ,ALLOCATABLE :: outside(:,:) + LOGICAL ,ALLOCATABLE :: outside_tmp(:,:) + LOGICAL ,POINTER :: out(:,:) +! REAL(rstd) :: pseudoalt(llm) + + INTEGER :: nbp_phys, nbp_phys_glo + +!$OMP PARALLEL + CALL allocate_field(f_bounds_lon,field_t,type_real,6) + CALL allocate_field(f_bounds_lat,field_t,type_real,6) + CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in") + CALL allocate_field(f_pks,field_t,type_real) + CALL allocate_field(f_pk,field_t,type_real,llm) + CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in") + CALL allocate_field(f_theta,field_t,type_real,llm) + CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in") + CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in") + CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in") + CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in") + CALL allocate_field(f_vort,field_z,type_real,llm,name="vort_in") + CALL allocate_field(f_vortc,field_t,type_real,llm,name="vortc_in") + CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out") + CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out") + CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out") + CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out") + CALL allocate_field(f_dps,field_t,type_real,name="dps_out") + CALL allocate_field(f_duc,field_t,type_real,3,llm) + + CALL init_message(f_dps,req_i0,req_dps0) + CALL init_message(f_dulon,req_i0,req_dulon0) + CALL init_message(f_dulat,req_i0,req_dulat0) + CALL init_message(f_dTemp,req_i0,req_dTemp0) + CALL init_message(f_dq,req_i0,req_dq0) +!$OMP END PARALLEL + + nbp_phys=0 + DO ind=1,ndomain + CALL swap_dimensions(ind) + + ALLOCATE(outside(ii_begin:ii_end,jj_begin:jj_end)) ! for limited area : don't take cells arround the border + ALLOCATE(outside_tmp(ii_begin-1:ii_end+1,jj_begin-1:jj_end+1)) ! for limited area : don't take cells arround the border + out=>domain(ind)%outside + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + outside(i,j)= out(i+1,j) .OR. & ! right + out(i,j+1 ) .OR. & ! rup + out(i-1 ,j+1) .OR. & ! lup + out(i-1 ,j) .OR. & !left + out(i ,j-1) .OR. & !ldown + out(i+1,j-1) !rdown + ENDDO + ENDDO + + outside_tmp=.FALSE. + outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside + + DO h=1,lam_halo_scheme-1 ! do not compute physic on limited area halo + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + outside(i,j) = outside_tmp(i,j) .OR. & + outside_tmp(i+1,j) .OR. & ! right + outside_tmp(i,j+1 ) .OR. & ! rup + outside_tmp(i-1 ,j+1) .OR. & ! lup + outside_tmp(i-1 ,j) .OR. & !left + outside_tmp(i ,j-1) .OR. & !ldown + outside_tmp(i+1,j-1) !rdown + ENDDO + ENDDO + outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside + ENDDO + + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + IF (domain(ind)%own(i,j) .AND. .NOT.outside(i,j)) nbp_phys=nbp_phys+1 + ENDDO + ENDDO + DEALLOCATE(outside) + DEALLOCATE(outside_tmp) + ENDDO + + +!initialize LMDZ5 physic mpi decomposition + CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr) + CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa) + + DO ind=1,ndomain + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + bounds_lon=f_bounds_lon(ind) + bounds_lat=f_bounds_lat(ind) + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1)) + CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2)) + CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3)) + CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4)) + CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5)) + CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6)) + ENDDO + ENDDO + ENDDO + + CALL init_infotrac_phy + + +!$OMP PARALLEL + CALL initialize_physics_omp + CALL param_gravity_wave +!$OMP END PARALLEL + + CALL xios_set_context + + END SUBROUTINE initialize_physics + + + SUBROUTINE initialize_physics_omp + USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz +! from dynamico + USE domain_mod + USE dimensions + USE mpi_mod + USE mpipara + USE disvert_mod + USE earth_const, ONLY: scale_height + USE xios_mod + USE time_mod , init_time_icosa=> init_time + USE omp_para + +! from LMDZ + USE mod_grid_phy_lmdz, ONLY : unstructured, klon_glo + USE mod_phys_lmdz_para, ONLY: klon_omp, reduce_min_lmdz => reduce_min , gather_lmdz => gather , bcast_lmdz => bcast + USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time + USE transfert_mod +! USE physics_distribution_mod, ONLY : init_physics_distribution + USE geometry_mod, ONLY : init_geometry + USE vertical_layers_mod, ONLY : init_vertical_layers + USE infotrac_phy, ONLY : init_infotrac_phy + USE inifis_mod, ONLY : inifis + USE readTracFiles_mod, ONLY: trac_type, isot_type + USE tracer_icosa_mod, ONLY : tracs + USE readTracFiles_mod, ONLY: delPhase +! USE phyaqua_mod, ONLY : iniaqua + + + IMPLICIT NONE + + + + INTEGER :: ind,i,j,k,ij,pos + REAL(rstd),POINTER :: bounds_lon(:,:) + REAL(rstd),POINTER :: bounds_lat(:,:) + + REAL(rstd),ALLOCATABLE :: latfi(:) + REAL(rstd),ALLOCATABLE :: lonfi(:) + REAL(rstd),ALLOCATABLE :: airefi(:) + REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:) + REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:) + REAL(rstd),ALLOCATABLE :: ind_cell_glo_r(:) + INTEGER, ALLOCATABLE :: ind_cell_glo(:) + INTEGER, ALLOCATABLE :: ind_cell_glo_tot(:) + INTEGER, ALLOCATABLE :: cell_glo_tot(:) + INTEGER :: ncell_glo_tot + + REAL(rstd) :: pseudoalt(llm) + REAL(rstd) :: aps(llm) + REAL(rstd) :: bps(llm) + REAL(rstd) :: scaleheight + + INTEGER :: run_length + REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!! + INTEGER :: annee_ref + INTEGER :: day_ref + INTEGER :: day_ini + REAL :: start_time + REAL :: physics_timestep + + ! Tracer stuff (SAVEd when needed to be OpenMP shared) + INTEGER :: nq + INTEGER :: nqo, nbtr, nbtr_inca + CHARACTER(len=256) :: type_trac + INTEGER,ALLOCATABLE :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it + INTEGER,ALLOCATABLE :: pbl_flg(:) ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it + CHARACTER(len=8),ALLOCATABLE :: solsym(:) ! tracer name from inca + + TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:) + + INTEGER :: iflag_phys + + INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema + INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema + + TYPE(trac_type) :: tracers_ico2lmdz(nqtot) !=== TRACERS DESCRIPTORS VECTOR + TYPE(isot_type) :: isotopes_ico2lmdz(1) !=== ISOTOPES PARAMETERS VECTOR + INTEGER :: iq + + CHARACTER(LEN=3) :: descrq(30) !--- Advection scheme description tags + logical, save :: first = .TRUE. + + + CALL init_distrib_icosa_lmdz + + ALLOCATE(latfi(klon_omp)) + ALLOCATE(lonfi(klon_omp)) + ALLOCATE(airefi(klon_omp)) + ALLOCATE(bounds_latfi(klon_omp,6)) + ALLOCATE(bounds_lonfi(klon_omp,6)) + ALLOCATE(ind_cell_glo_r(klon_omp)) + ALLOCATE(ind_cell_glo(klon_omp)) + ALLOCATE(q_ave(llm,nqtot)) + + CALL transfer_icosa_to_lmdz(geom%lat_i,latfi) + CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi) + CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi) + CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi) + CALL transfer_icosa_to_lmdz(geom%Ai,airefi) + + CALL allocate_field(f_ind_cell_glo,field_t,type_real) + + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j) + ENDDO + ENDDO + ENDDO + + + CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo_r) + CALL deallocate_field(f_ind_cell_glo) + ind_cell_glo=INT(ind_cell_glo_r) + DEALLOCATE(ind_cell_glo_r) + + + CALL reduce_min_lmdz(MINVAL(-ind_cell_glo),ncell_glo_tot) ! reduce_max does not exist in lmdz, use reduce_min + CALL bcast_lmdz(ncell_glo_tot) + ncell_glo_tot=-ncell_glo_tot + ALLOCATE(cell_glo_tot(0:ncell_glo_tot)) + ALLOCATE(ind_cell_glo_tot(klon_glo)) + CALL gather_lmdz(ind_cell_glo,ind_cell_glo_tot) + CALL bcast_lmdz(ind_cell_glo_tot) + + cell_glo_tot=-1 + DO i=1,klon_glo + cell_glo_tot(ind_cell_glo_tot(i))= 0 + ENDDO + + pos=0 + DO i=0,ncell_glo_tot + IF (cell_glo_tot(i)/=-1) THEN + cell_glo_tot(i) = pos + pos=pos + 1 + ENDIF + ENDDO + + DO i=1,klon_omp + ind_cell_glo(i)=cell_glo_tot(ind_cell_glo(i)) + ENDDO + + ind_cell_glo = ind_cell_glo + 1 ! lmdz expect global indices begining to 1 not 0 + + +! CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_mpi,1,MPI_INTEGER,comm_icosa,ierr) +! +! displ(0)=0 +! DO i=1,mpi_size-1 +! displ(i)=displ(i-1)+ncell_mpi(i-1) +! ENDDO +! +! ALLOCATE(ind_glo_tot(ncell_tot)) +! ALLOCATE(cell_glo_tot(0:ncell_glo_tot-1)) +! +! cell_glo_tot(:)= -1 +! CALL MPI_ALLGATHERV(ind_glo, ncell, MPI_INTEGER, ind_glo_tot, ncell_mpi, displ, MPI_INTEGER, comm_icosa,ierr) +! +! DO i=1,ncell_tot +! cell_glo_tot(ind_glo_tot(i))= 0 +! ENDDO +! +! ncell_glo=0 +! DO i=0,ncell_glo_tot-1 +! IF (cell_glo_tot(i)/=-1) THEN +! cell_glo_tot(i) = ncell_glo +! ncell_glo=ncell_glo + 1 +! ENDIF +! ENDDO +! +! DO i=1,ncell +! ind_glo(i)=cell_glo_tot(ind_glo(i)) +! ENDDO + + + CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, ind_cell_glo) + + scaleheight=scale_height/1000. ! Atmospheric scale height (km) + aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1)) + bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1)) + pseudoalt(:)=-scaleheight*log(presnivs(:)/preff) + CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,presinter,pseudoalt) + + ! Initialize tracer names, numbers, etc. for physics + !Config Key = type_trac + !Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS + !Config Def = lmdz + !Config Help = + !Config 'lmdz' = pas de couplage, pur LMDZ + !Config 'lmdz|inca' = model de chime INCA + !Config 'lmdz|repr' = model de chime REPROBUS + type_trac = 'lmdz' + CALL getin('type_trac',type_trac) + + descrq( 1: 2) = ['LMV','BAK'] + descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH',' ','PPM','PPS','PPP',' ','SLP'] + descrq(30) = 'PRA' + + nqo = 0 + DO iq=1,nqtot + + tracers_ico2lmdz(iq)%name = tracs(iq)%name + + tracers_ico2lmdz(iq)%gen0Name = tracs(iq)%name + tracers_ico2lmdz(iq)%phase = tracs(iq)%phase + tracers_ico2lmdz(iq)%iadv = tracs(iq)%iadv + + IF (tracs(iq)%component .eq. "dynamico") then + tracers_ico2lmdz(iq)%component='lmdz' + ELSE + tracers_ico2lmdz(iq)%component=tracs(iq)%component + ENDIF + + if (tracers_ico2lmdz(iq)%iadv .ne. 0 ) tracers_ico2lmdz(iq)%isAdvected=.true. + + tracers_ico2lmdz(iq)%longName = tracers_ico2lmdz(iq)%name + IF(tracers_ico2lmdz(iq)%iadv > 0) tracers_ico2lmdz(iq)%longName=TRIM(tracers_ico2lmdz(iq)%name)//descrq(tracers_ico2lmdz(iq)%iadv) + + tracers_ico2lmdz(iq)%isInPhysics= delPhase(tracers_ico2lmdz(iq)%gen0Name) /= 'H2O' .OR. tracers_ico2lmdz(iq)%component /= 'lmdz' + tracers_ico2lmdz(iq)%iGeneration = 0 + + ENDDO + + nqo = COUNT(delPhase(tracs(:)%name) == 'H2O' .AND. tracers_ico2lmdz(:)%component == 'lmdz') !--- Number of water phases + + isotopes_ico2lmdz(1)%parent='H2O' + isotopes_ico2lmdz(1)%phase='gls' + isotopes_ico2lmdz(1)%nphas=3 + + + nbtr=nqtot-nqo + + ALLOCATE(conv_flg(nbtr)) + ALLOCATE(pbl_flg(nbtr)) + + conv_flg(:) = 1 ! convection activated for all tracers + pbl_flg(:) = 1 ! boundary layer activated for all tracers + + + CALL init_infotrac_phy + + ! Initialize physical constant + day_length=86400 + CALL getin('day_length',day_length) + CALL inifis(day_length,radius,g,kappa*cpp,cpp) + + + + ! init time + annee_ref=2015 + CALL getin("anneeref",annee_ref) + + day_ref=1 + CALL getin("dayref",day_ref) + + physics_timestep=dt*itau_physics + run_length=itaumax*dt + ndays=NINT(run_length/day_length) + + day_ini=INT(itau0*dt/day_length)+day_ref + start_time= itau0*dt/day_length-INT(itau0*dt/day_length) + + CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, int(ndays), physics_timestep) + +! Additional initializations for aquaplanets +! CALL getin("iflag_phys",iflag_phys) +! IF (iflag_phys>=100) THEN +! CALL iniaqua(klon_omp, iflag_phys) +! END IF + + + +#ifdef INCA + CONTAINS + + SUBROUTINE init_chem_trac() + IMPLICIT NONE + + CALL Init_chem_inca_trac(nbtr) + + END SUBROUTINE init_chem_trac + + SUBROUTINE init_chem_transport() + + IMPLICIT NONE + + CALL init_transport(solsym, conv_flg,pbl_flg, hadv_inca, vadv_inca) + + END SUBROUTINE init_chem_transport + + +#else + CONTAINS + SUBROUTINE init_chem_trac() + IMPLICIT NONE + + END SUBROUTINE init_chem_trac + + SUBROUTINE init_chem_transport() + + IMPLICIT NONE + + END SUBROUTINE init_chem_transport + + +#endif + + + + END SUBROUTINE initialize_physics_omp + + + + + SUBROUTINE physics + USE icosa + USE time_mod + USE disvert_mod + USE transfert_mod + USE mpipara + USE omp_para + USE xios_mod + USE wxios + USE trace + USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa + USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q + USE write_field_mod + USE checksum_mod + USE vorticity_mod + USE tracer_icosa_mod + USE compute_transport_mod + USE, INTRINSIC :: IEEE_ARITHMETIC + +! from LMDZ + USE mod_phys_lmdz_omp_data, ONLY: klon_omp + USE geometry_mod, ONLY : cell_area + USE physiq_mod, ONLY: physiq + USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave + USE mod_phys_lmdz_para, ONLY: reduce_sum, reduce_min + USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ + USE write_field_phy + IMPLICIT NONE + + REAL(rstd),POINTER :: phis(:) + REAL(rstd),POINTER :: ps(:) + REAL(rstd),POINTER :: theta_rhodz(:,:,:) + REAL(rstd),POINTER :: u(:,:) + REAL(rstd),POINTER :: wflux(:,:) + REAL(rstd),POINTER :: q(:,:,:) + REAL(rstd),POINTER :: p(:,:) + REAL(rstd),POINTER :: pks(:) + REAL(rstd),POINTER :: pk(:,:) + REAL(rstd),POINTER :: p_layer(:,:) + REAL(rstd),POINTER :: theta(:,:) + REAL(rstd),POINTER :: phi(:,:) + REAL(rstd),POINTER :: Temp(:,:) + REAL(rstd),POINTER :: ulon(:,:) + REAL(rstd),POINTER :: ulat(:,:) + REAL(rstd),POINTER :: vort(:,:) + REAL(rstd),POINTER :: vortc(:,:) + REAL(rstd),POINTER :: dulon(:,:) + REAL(rstd),POINTER :: dulat(:,:) + REAL(rstd),POINTER :: dTemp(:,:) + REAL(rstd),POINTER :: dq(:,:,:) + REAL(rstd),POINTER :: dps(:) + REAL(rstd),POINTER :: duc(:,:,:) + + + INTEGER :: ind,l,k, iq, iq_child + + REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:) +!$OMP THREADPRIVATE(ps_phy) + REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:) +!$OMP THREADPRIVATE(p_phy) + REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:) +!$OMP THREADPRIVATE(p_layer_phy) + REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:) +!$OMP THREADPRIVATE(Temp_phy) + REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:) +!$OMP THREADPRIVATE(phis_phy) + REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:) +!$OMP THREADPRIVATE(phi_phy) + REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:) +!$OMP THREADPRIVATE(ulon_phy) + REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:) +!$OMP THREADPRIVATE(ulat_phy) + REAL(rstd),ALLOCATABLE,SAVE :: rot_phy(:,:) +!$OMP THREADPRIVATE(rot_phy) + REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:) +!$OMP THREADPRIVATE(q_phy) + REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:) +!$OMP THREADPRIVATE(wflux_phy) + REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:) +!$OMP THREADPRIVATE(dulon_phy) + REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:) +!$OMP THREADPRIVATE(dulat_phy) + REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:) +!$OMP THREADPRIVATE(dTemp_phy) + REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:) +!$OMP THREADPRIVATE(dq_phy) + REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:) +!$OMP THREADPRIVATE(dps_phy) + REAL(rstd) :: dtphy + LOGICAL :: debut + LOGICAL :: lafin + LOGICAL,SAVE :: first=.TRUE. +!$OMP THREADPRIVATE(first) + INTEGER, SAVE :: iter +!$OMP THREADPRIVATE(iter) + REAL(rstd),ALLOCATABLE,SAVE :: q_phy_old(:,:,:) + REAL :: qmin_omp,qmin_glo + REAL :: q_ave_bis(llm,nqtot) + + IF(first) THEN + debut=.TRUE. + iter = 1 + ELSE + debut=.FALSE. + iter=iter+1 + ENDIF + + + IF(it-itau0>=itaumax) THEN + lafin=.TRUE. + ELSE + lafin=.FALSE. + ENDIF + + IF (first) THEN + first=.FALSE. + CALL init_message(f_u,req_e1_vect,req_u) + CALL init_message(f_vort,req_z1_scal,req_z) + ALLOCATE(ps_phy(klon_omp)) + ALLOCATE(p_phy(klon_omp,llm+1)) + ALLOCATE(p_layer_phy(klon_omp,llm)) + ALLOCATE(Temp_phy(klon_omp,llm)) + ALLOCATE(phis_phy(klon_omp)) + ALLOCATE(phi_phy(klon_omp,llm)) + ALLOCATE(ulon_phy(klon_omp,llm)) + ALLOCATE(ulat_phy(klon_omp,llm)) + ALLOCATE(rot_phy(klon_omp,llm)) + ALLOCATE(q_phy(klon_omp,llm,nqtot)) + ALLOCATE(wflux_phy(klon_omp,llm)) + ALLOCATE(dulon_phy(klon_omp,llm)) + ALLOCATE(dulat_phy(klon_omp,llm)) + ALLOCATE(dTemp_phy(klon_omp,llm)) + ALLOCATE(dq_phy(klon_omp,llm,nqtot)) + ALLOCATE(dps_phy(klon_omp)) + ALLOCATE(q_phy_old(klon_omp,llm,nqtot)) +!$OMP BARRIER + +! CALL param_gravity_wave + ENDIF + + +!$OMP MASTER +! CALL update_calendar(it) +!$OMP END MASTER +!$OMP BARRIER + dtphy=itau_physics*dt + + + + CALL transfert_message(f_u,req_u) + DO ind=1,ndomain + IF (assigned_domain(ind)) THEN + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + u=f_u(ind) + vort=f_vort(ind) + CALL compute_vorticity(u,vort) + ENDIF + ENDDO + + CALL transfert_message(f_vort,req_z) + + + DO ind=1,ndomain + CALL swap_dimensions(ind) + IF (assigned_domain(ind)) THEN + CALL swap_geometry(ind) + + phis=f_phis(ind) + ps=f_ps(ind) + theta_rhodz=f_theta_rhodz(ind) + u=f_u(ind) + q=f_q(ind) + wflux=f_wflux(ind) + p=f_p(ind) + pks=f_pks(ind) + pk=f_pk(ind) + p_layer=f_p_layer(ind) + theta=f_theta(ind) + phi=f_phi(ind) + Temp=f_Temp(ind) + ulon=f_ulon(ind) + ulat=f_ulat(ind) + vort=f_vort(ind) + vortc=f_vortc(ind) + + CALL grid_icosa_to_physics + + ENDIF + ENDDO + +!$OMP BARRIER +!$OMP MASTER + CALL SYSTEM_CLOCK(start_clock) +!$OMP END MASTER + CALL trace_start("physic") +! CALL trace_off() + + +! CALL writeField("p_in",f_p) +! CALL writeField("p_layer_in",f_p_layer) +! CALL writeField("phi_in",f_phi) +! CALL writeField("phis_in",f_phis) +! CALL writeField("ulon_in",f_ulon) +! CALL writeField("ulat_in",f_ulat) +! CALL writeField("Temp_in",f_Temp) +! CALL writeField("q_in",f_q) +! CALL writeField("wflux_in",f_wflux) +! CALL writeField("vortc",f_vortc) + +! CALL checksum(f_p) +! CALL checksum(f_p_layer) +! CALL checksum(f_phi) +! CALL checksum(f_phis) +! CALL checksum(f_ulon) +! CALL checksum(f_ulat) +! CALL checksum(f_Temp) +! CALL checksum(f_q) +! CALL checksum(f_wflux) + + CALL transfer_icosa_to_lmdz(f_p , p_phy) + CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy) + CALL transfer_icosa_to_lmdz(f_phi , phi_phy) + CALL transfer_icosa_to_lmdz(f_phis , phis_phy ) + CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy ) + CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy) + CALL transfer_icosa_to_lmdz(f_vortc , rot_phy) + CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy) + CALL transfer_icosa_to_lmdz(f_q , q_phy) + CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy) + + DO l=1,llm + wflux_phy(:,l) = - wflux_phy(:,l)*cell_area(:) + phi_phy(:,l)=phi_phy(:,l)-phis_phy(:) + ENDDO + + DO iq=1, nqtot + IF (.NOT. tracers(iq)%has_parent) THEN + DO iq_child=1,tracers(iq)%nb_children + CALL tracers_dyn2phys(q_phy, tracers(iq)%children(iq_child), iq) + ENDDO + ENDIF + + ENDDO + + CALL wxios_set_context() + + ! Ehouarn: rot_phy() not implemented!! Set it to zero for now +! rot_phy(:,:)=0 + CALL physiq(klon_omp, llm, debut, lafin, dtphy, & + p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, & + ulon_phy, ulat_phy, rot_phy, Temp_phy, q_phy, wflux_phy, & + dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy) + + CALL compute_q_ave(q_phy, q_ave) + + DO iq=1, nqtot + q_phy(:,:,iq) = q_phy(:,:,iq) + dtphy*dq_phy(:,:,iq) + ENDDO + +! IF (is_master) PRINT*, "==================================================" +! IF (is_master) PRINT*, "==================================================" +! CALL compute_q_ave(q_phy, q_ave_bis) +! DO iq=1, nqtot +! qmin_omp = MINVAL(q_phy(:,:,iq)) +! CALL reduce_min(qmin_omp, qmin_glo) +! IF (is_master) PRINT*, "after 0 : MINVAL of ",TRIM(tracers(iq)%name), qmin_glo +! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave_bis(:,iq) +! ENDDO + + CALL compute_qminimum(q_phy, q_ave, p_phy) + +! DO iq=1, nqtot +! IF (.NOT. tracers(iq)%has_parent) THEN +! DO iq_child=1,tracers(iq)%nb_children +! CALL tracers_normalize(q_phy, tracers(iq)%children(iq_child), iq) +! ENDDO +! DO l=1,llm +! DO k=1,klon_omp +! IF ( q_phy(k,l,iq) < q_epsilon*q_ave(l,iq)) q_phy(k,l,iq) = q_epsilon*q_ave(l,iq) +! ENDDO +! ENDDO +! ENDIF +! ENDDO + + +! IF (is_master) PRINT*, "==================================================" +! IF (is_master) PRINT*, "==================================================" +! +! CALL compute_q_ave(q_phy, q_ave) +! +! DO iq=1, nqtot +! qmin_omp = MINVAL(q_phy(:,:,iq)) +! CALL reduce_min(qmin_omp, qmin_glo) +! IF (is_master) PRINT*, "after 1 MINVAL of ",TRIM(tracers(iq)%name), qmin_glo +! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave(:,iq) +! IF (ANY(IEEE_IS_NAN(q_ave(:,iq)))) THEN +! CALL WriteField_phy("debug_"//TRIM(tracers(iq)%name),q_phy(:,:,iq),llm) +! IF (is_master) PRINT *,'NAN DETECTED' +! STOP +! ENDIF +! ENDDO + + DO iq=1, nqtot + IF (.NOT. tracers(iq)%has_parent) THEN + DO iq_child=1,tracers(iq)%nb_children + CALL tracers_phys2dyn(q_phy, tracers(iq)%children(iq_child), iq) + ENDDO + ENDIF + + ENDDO + + + CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon ) + CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat ) + CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp ) + CALL transfer_lmdz_to_icosa(q_phy , f_q ) + CALL transfer_lmdz_to_icosa(dps_phy , f_dps ) + +! CALL writeField("dulon_out",f_dulon) +! CALL writeField("dulat_out",f_dulat) +! CALL writeField("dTemp_out",f_dTemp) +! CALL writeField("dq_out",f_dq) +! CALL writeField("dps_out",f_dps) + +! CALL checksum(f_dulon) +! CALL checksum(f_dulat) +! CALL checksum(f_dTemp) +! CALL checksum(f_dq) +! CALL checksum(f_dps) + + CALL send_message(f_dps,req_dps0) + CALL send_message(f_dulon,req_dulon0) + CALL send_message(f_dulat,req_dulat0) + CALL send_message(f_dTemp,req_dTemp0) + CALL send_message(f_q,req_dq0) + + CALL wait_message(req_dps0) + CALL wait_message(req_dulon0) + CALL wait_message(req_dulat0) + CALL wait_message(req_dTemp0) + CALL wait_message(req_dq0) + + +! CALL trace_on() + CALL trace_end("physic") +!$OMP MASTER + CALL SYSTEM_CLOCK(stop_clock) + count_clock=count_clock+stop_clock-start_clock +!$OMP END MASTER + +!$OMP BARRIER + + DO ind=1,ndomain + CALL swap_dimensions(ind) + IF (assigned_domain(ind)) THEN + CALL swap_geometry(ind) + + theta_rhodz=f_theta_rhodz(ind) + u=f_u(ind) + q=f_q(ind) + ps=f_ps(ind) + dulon=f_dulon(ind) + dulat=f_dulat(ind) + Temp=f_temp(ind) + dTemp=f_dTemp(ind) + dq=f_dq(ind) + dps=f_dps(ind) + duc=f_duc(ind) + p=f_p(ind) + pks=f_pks(ind) + pk=f_pk(ind) + + CALL grid_physics_to_icosa + ENDIF + ENDDO + +!$OMP BARRIER + CALL xios_set_context + + + CONTAINS + + SUBROUTINE grid_icosa_to_physics + USE pression_mod + USE exner_mod + USE theta2theta_rhodz_mod + USE geopotential_mod + USE wind_from_lonlat_mod + USE omp_para + IMPLICIT NONE + + REAL(rstd) :: uc(3) + INTEGER :: i,j,ij,l + +! compute pression + + DO l = ll_begin,ll_endp1 + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + p(ij,l) = ap(l) + bp(l) * ps(ij) + ENDDO + ENDDO + ENDDO + +!$OMP BARRIER + +! compute exner + + IF (is_omp_first_level) THEN + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + pks(ij) = cpp * ( ps(ij)/preff ) ** kappa + ENDDO + ENDDO + ENDIF + + ! 3D : pk + DO l = ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa + ENDDO + ENDDO + ENDDO + +!$OMP BARRIER + +! compute theta, temperature and pression at layer + DO l = ll_begin, ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g) + Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp + p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) + ENDDO + ENDDO + ENDDO + + +!!! Compute geopotential + + ! for first layer + IF (is_omp_first_level) THEN + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) + ENDDO + ENDDO + ENDIF +!!-> implicit flush on phi(:,1) + +!$OMP BARRIER + + ! for other layers + DO l = ll_beginp1, ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + phi(ij,l) = 0.5 * ( theta(ij,l) + theta(ij,l-1) ) & + * ( pk(ij,l-1) - pk(ij,l) ) + ENDDO + ENDDO + ENDDO + +!$OMP BARRIER + + + IF (is_omp_first_level) THEN + DO l = 2, llm + DO j=jj_begin,jj_end +! ---> Bug compilo intel ici en openmp +! ---> Couper la boucle + IF (j==jj_end+1) PRINT*,"this message must not be printed" + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + phi(ij,l) = phi(ij,l)+ phi(ij,l-1) + ENDDO + ENDDO + ENDDO +! --> IMPLICIT FLUSH on phi --> non + ENDIF + +! compute wind centered lon lat compound + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + uc(:)=1/Ai(ij)* & + ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:)) & + + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:)) & + + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:)) & + + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:)) & + + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))& + + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:))) + ulon(ij,l)=sum(uc(:)*elon_i(ij,:)) + ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) + ENDDO + ENDDO + ENDDO + + +! compute centered vorticity + + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + vortc(ij,l) = Riv(ij,vup) * vort(ij+z_up,l) + & + Riv(ij,vlup) * vort(ij+z_lup,l) + & + Riv(ij,vldown)* vort(ij+z_ldown,l) + & + Riv(ij,vdown) * vort(ij+z_down,l) + & + Riv(ij,vrdown)* vort(ij+z_rdown,l) + & + Riv(ij,vrup) * vort(ij+z_rup,l) + ENDDO + ENDDO + ENDDO + + + +!$OMP BARRIER + END SUBROUTINE grid_icosa_to_physics + + + SUBROUTINE grid_physics_to_icosa + USE theta2theta_rhodz_mod + USE omp_para + IMPLICIT NONE + INTEGER :: i,j,ij,l,iq + + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:) + ENDDO + ENDDO + ENDDO + + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) + u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) + u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) + ENDDO + ENDDO + ENDDO + + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l) + ENDDO + ENDDO + ENDDO + +! DO iq=1,nqtot +! DO l=ll_begin,ll_end +! DO j=jj_begin,jj_end +! DO i=ii_begin,ii_end +! ij=(j-1)*iim+i +! q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq) +! ENDDO +! ENDDO +! ENDDO +! ENDDO + +!$OMP BARRIER + + IF (is_omp_first_level) THEN + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + ps(ij)=ps(ij)+ dtphy * dps(ij) + ENDDO + ENDDO + ENDIF + +! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0) + +! compute pression +!$OMP BARRIER + DO l = ll_begin,ll_endp1 + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + p(ij,l) = ap(l) + bp(l) * ps(ij) + ENDDO + ENDDO + ENDDO + +!$OMP BARRIER + +! compute exner + + IF (is_omp_first_level) THEN + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + pks(ij) = cpp * ( ps(ij)/preff ) ** kappa + ENDDO + ENDDO + ENDIF + + ! 3D : pk + DO l = ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa + ENDDO + ENDDO + ENDDO + +!$OMP BARRIER + +! compute theta, temperature and pression at layer + DO l = ll_begin, ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp ) + ENDDO + ENDDO + ENDDO + + END SUBROUTINE grid_physics_to_icosa + + RECURSIVE SUBROUTINE tracers_dyn2phys(q, iq, iq_parent) + USE tracer_icosa_mod + IMPLICIT NONE + REAL(rstd),INTENT(INOUT) :: q(:,:,:) + INTEGER,INTENT(IN) :: iq, iq_parent + INTEGER :: iq_child + + q(:,:,iq) = q(:,:,iq)*q(:,:,iq_parent) + + DO iq_child=1,tracers(iq)%nb_children + CALL tracers_dyn2phys(q, tracers(iq)%children(iq_child), iq) + ENDDO + + END SUBROUTINE tracers_dyn2phys + + RECURSIVE SUBROUTINE tracers_phys2dyn(q, iq, iq_parent) + USE tracer_icosa_mod + USE compute_transport_mod + USE mod_phys_lmdz_para + USE mod_phys_lmdz_omp_data, ONLY: klon_omp + USE dimphy + USE mod_grid_phy_lmdz + IMPLICIT NONE + REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot) + INTEGER,INTENT(IN) :: iq, iq_parent + INTEGER :: iq_child + INTEGER :: lij + + DO iq_child=1,tracers(iq)%nb_children + CALL tracers_phys2dyn(q, tracers(iq)%children(iq_child), iq) + ENDDO + + q(:,:,iq) = q(:,:,iq) / q(:,:,iq_parent) + + END SUBROUTINE tracers_phys2dyn + + RECURSIVE SUBROUTINE tracers_normalize(q, iq, iq_parent) + USE tracer_icosa_mod + USE compute_transport_mod + USE mod_phys_lmdz_para + USE mod_phys_lmdz_omp_data, ONLY: klon_omp + USE dimphy + USE mod_grid_phy_lmdz + IMPLICIT NONE + REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot) + INTEGER,INTENT(IN) :: iq, iq_parent + INTEGER :: iq_child + INTEGER :: l, k + + DO iq_child=1,tracers(iq)%nb_children + CALL tracers_normalize(q, tracers(iq)%children(iq_child), iq) + ENDDO + + DO l=1,klev + DO k=1,klon + IF ( q(k,l,iq_parent) >= q_epsilon*q_ave(l, iq_parent) ) THEN + IF ( q(k,l,iq) < q_epsilon*q_ave(l, iq_parent) ) THEN + q(k,l,iq) = q_ave(l, iq)/q_ave(l, iq_parent) * q(k,l,iq_parent) + ENDIF + ELSE + q(k,l,iq) = q_ave(l, iq) * q_epsilon + ENDIF + ENDDO + ENDDO + + END SUBROUTINE tracers_normalize + + + SUBROUTINE compute_q_ave(q, q_ave) + USE dimphy + USE tracer_icosa_mod + USE mod_phys_lmdz_para + USE mod_grid_phy_lmdz + IMPLICIT NONE + REAL(rstd),INTENT(IN) :: q(klon, klev ,nqtot) + REAL(rstd),INTENT(OUT) :: q_ave(klev, nqtot) + REAL(rstd) :: q_ave_omp(klev, nqtot) + INTEGER :: iq, l , k + + DO iq=1,nqtot + DO l=1,klev + q_ave_omp(l,iq) = 0. + DO k=1,klon + q_ave_omp(l,iq) = q_ave_omp(l,iq) + q(k,l,iq) + ENDDO + ENDDO + ENDDO + CALL reduce_sum(q_ave_omp, q_ave) + CALL bcast(q_ave) + q_ave(:,:) = q_ave(:,:) / klon_glo + WHERE (q_ave < q_minimum) q_ave = q_minimum + + END SUBROUTINE compute_q_ave + + SUBROUTINE compute_qminimum(q, q_ave, p ) + USE dimphy + USE tracer_icosa_mod + IMPLICIT NONE + REAL(rstd),INTENT(INOUT) :: q(klon, klev ,nqtot) + REAL(rstd),INTENT(IN) :: q_ave(klev, nqtot) + REAL(rstd),INTENT(IN) :: p(klon, klev) + REAL(rstd) :: deltap(klon, klev) + INTEGER :: iq, iq_g + INTEGER :: iq_H2O_g, iq_H2O_l, iq_H216O_g, iq_H216O_l , iq_H218O_g, iq_H218O_l, iq_HDO_g, iq_HDO_l + + DO iq=1, nqtot + SELECT CASE (TRIM(tracers(iq)%name)) + CASE ("H2O_g") ; iq_H2O_g = iq + CASE ("H2O_l") ; iq_H2O_l = iq + CASE ("H216O_g") ; iq_H216O_g = iq + CASE ("H216O_l") ; iq_H216O_l = iq + CASE ("H218O_g") ; iq_H218O_g = iq + CASE ("H218O_l") ; iq_H218O_l = iq + CASE ("HDO_g") ; iq_HDO_g = iq + CASE ("HDO_l") ; iq_HDO_l = iq + CASE DEFAULT ; CYCLE + END SELECT + ENDDO + + DO iq=1, nqtot + IF (iq==iq_H2O_l) THEN + iq_g = iq_H2O_g + ELSE IF (iq==iq_H216O_l) THEN + iq_g = iq_H216O_g + ELSE IF (iq==iq_H218O_l) THEN + iq_g = iq_H218O_g + ELSE IF (iq==iq_HDO_l) THEN + iq_g = iq_HDO_g + ELSE + CYCLE + ENDIF + + DO l=1,klev + DO k=1,klon + IF ( q(k,l,iq) <= q_epsilon*q_ave(l,iq)) THEN + q(k,l,iq_g) = q(k,l,iq_g) - (q_epsilon*q_ave(l,iq) - q(k,l,iq)) + q(k,l,iq) = q_epsilon*q_ave(l,iq) + ENDIF + ENDDO + ENDDO + ENDDO + + DO l=1,klev + DO k=1,klon + deltap(k,l) = p(k,l)-p(k,l+1) + ENDDO + ENDDO + + DO iq=1, nqtot + DO l=klev,2 + DO k=1,klon + IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN + q(k,l-1,iq) = q(k,l-1,iq) - (q_epsilon*q_ave(l,iq) - q(k,l,iq))*deltap(k,l)/deltap(k,l-1) + q(k,l,iq) = q_epsilon*q_ave(l,iq) + ENDIF + ENDDO + ENDDO + ENDDO + + DO iq=1, nqtot + DO l=1,klev-1 + DO k=1,klon + IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN + q(k,l+1,iq) = q(k,l+1,iq) - (q_epsilon*q_ave(l,iq) - q(k,l,iq))*deltap(k,l)/deltap(k,l+1) + q(k,l,iq) = q_epsilon*q_ave(l,iq) + ENDIF + ENDDO + + DO k=1,klon + IF (q(k,klev,iq) < q_epsilon*q_ave(klev,iq)) q(k,klev,iq) = q_epsilon*q_ave(klev,iq) + ENDDO + ENDDO + ENDDO + + END SUBROUTINE compute_qminimum + + + END SUBROUTINE physics + +END MODULE interface_icosa_lmdz_mod