Ignore:
Timestamp:
Nov 10, 2023, 4:48:52 PM (14 months ago)
Author:
jbclement
Message:

PEM:
Correction of the reading of the PCM data (it did not work if no slope was used) + some minor related cleanings.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3114 r3122  
    130130== 03/11/2023 == JBC
    131131Following r3113, addition of 'nqsoil' and 'qsoil' in the arguments of the subroutines 'phyetat0' and 'physdem1' to be able to compile.
     132
     133== 10/11/2023 == JBC
     134Correction of the reading of the PCM data (it did not work if no slope was used) + some minor related cleanings.

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

    r3114 r3122  
    604604!------------------------
    605605call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
    606               porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,          &
    607               tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                      &
    608               qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,             &
     606              porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,       &
     607              tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                   &
     608              qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,          &
    609609              co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
    610610
    … …  
    717717        do t = 1,timelen
    718718            q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
    719                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
    720                                 + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
    721                                 -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
     719                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
     720                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
     721                                -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/                     &
    722722                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
    723723            if (q_co2_PEM_phys(ig,t) < 0) then
    … …  
    861861
    862862        deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
     863
     864! II_d.3 Update the ice table
    863865        write(*,*) "Compute ice table"
    864 
    865 ! II_d.3 Update the ice table
    866866        porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
    867867        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
    868 
    869868        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    870869
     870! II_d.4 Update the soil thermal properties
    871871        write(*,*) "Update soil propreties"
    872 
    873 ! II_d.4 Update the soil thermal properties
    874872        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
    875873
    … …  
    932930    i_myear = i_myear + dt_pem
    933931
    934     write(*,*) "Checking all the stopping criterion."
     932    write(*,*) "Checking all the stopping criteria..."
    935933    if (STOPPING_water) then
    936934        write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
    … …  
    11361134             float(day_ini),0.,nslope,def_slope,subslope_dist)
    11371135call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
    1138              TI_PEM, porefillingice_depth,porefillingice_thickness,         &
     1136             TI_PEM, porefillingice_depth,porefillingice_thickness,      &
    11391137             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    11401138write(*,*) "restartpem.nc has been written"

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

    r3096 r3122  
    77implicit none
    88
    9 character(256) :: msg, var, modname ! for reading
    10 integer        :: fID, vID          ! for reading
     9character(13), parameter :: modname = 'read_data_PCM'
     10character(256)           :: msg      ! for reading
     11integer                  :: fID, vID ! for reading
    1112
    1213!=======================================================================
    … …  
    1415!=======================================================================
    1516
    16 SUBROUTINE read_data_PCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
     17SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,          &
    1718                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
    1819                         watersurf_density_ave,watersoil_density)
    … …  
    3536! Arguments:
    3637
    37 character(len = *), intent(in) :: fichnom                             !--- FILE NAME
     38character(len = *), intent(in) :: filename                            ! File name
    3839integer,            intent(in) :: timelen                             ! number of times stored in the file
    3940integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    … …  
    5354real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
    5455real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
    55 !===============================================================================
    56 !   Local Variables
     56!=======================================================================
     57! Local Variables
    5758character(12)                   :: start_file_type = "earth" ! default start file type
    5859real, dimension(:), allocatable :: time      ! times stored in start
    … …  
    6061integer                         :: edges(4), corner(4)
    6162integer                         :: i, j, l, t          ! loop variables
    62 real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
     63real                            :: A, B, mmean         ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
    6364integer                         :: islope              ! loop for variables
    6465character(2)                    :: num                 ! for reading sloped variables
    6566real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
    66 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
    67 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
     67real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
     68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer [mol/m^3]
    6869real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
    6970real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
    … …  
    7980real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
    8081real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
    81 
    8282!-----------------------------------------------------------------------
    83 modname="read_data_PCM"
    84 
    85 A = (1/m_co2 - 1/m_noco2)
    86 B = 1/m_noco2
    87 
    88 write(*,*) "Opening ", fichnom, "..."
    89 
    90 !  Open initial state NetCDF file
    91 var = fichnom
    92 call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    93 
    94 write(*,*) "Downloading data for vmr co2..."
    95 
    96 call get_var3("co2_layer1"   ,q_co2_dyn)
    97 
    98 write(*,*) "Downloading data for vmr co2 done"
    99 write(*,*) "Downloading data for vmr h20..."
    100 
    101 call get_var3("h2o_layer1"   ,q_h2o_dyn)
    102 
    103 write(*,*) "Downloading data for vmr h2o done"
    104 write(*,*) "Downloading data for surface pressure ..."
    105 
    106 call get_var3("ps"   ,ps_GCM)
    107 
    108 write(*,*) "Downloading data for surface pressure done"
    109 write(*,*) "nslope=", nslope
    110 
    111 if (nslope > 1) then
    112     write(*,*) "Downloading data for co2ice_slope ..."
    113 
     83! Open the NetCDF file
     84write(*,*) "Opening "//filename//"..."
     85call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
     86
     87! Dowload the data from the file
     88call get_var3("ps",ps_GCM)
     89write(*,*) "Data for surface pressure downloaded!"
     90
     91call get_var3("co2_layer1",q_co2_dyn)
     92write(*,*) "Data for vmr co2 downloaded!"
     93
     94call get_var3("h2o_layer1",q_h2o_dyn)
     95write(*,*) "Data for vmr h2o downloaded!"
     96
     97if (nslope == 1) then ! There is no slope
     98    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
     99    write(*,*) "Data for co2_ice downloaded!"
     100
     101    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
     102    write(*,*) "Data for h2o_ice_s downloaded!"
     103
     104#ifndef CPP_STD
     105    call get_var3("watercap",watercap(:,:,1,:))
     106    write(*,*) "Data for watercap downloaded!"
     107#endif
     108
     109    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
     110    write(*,*) "Data for tsurf downloaded!"
     111
     112#ifndef CPP_STD
     113    if (soil_pem) then
     114        call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
     115        write(*,*) "Data for soiltemp downloaded!"
     116
     117        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
     118        write(*,*) "Data for waterdensity_soil downloaded!"
     119
     120        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,islope,:))
     121        write(*,*) "Data for waterdensity_surface downloaded!"
     122    endif !soil_pem
     123#endif
     124else ! We use slopes
    114125    do islope = 1,nslope
    115126        write(num,'(i2.2)') islope
    116127        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    117128    enddo
    118    
    119     write(*,*) "Downloading data for co2ice_slope done"
    120     write(*,*) "Downloading data for h2o_ice_s_slope ..."
    121    
     129    write(*,*) "Data for co2_ice downloaded!"
     130
    122131    do islope = 1,nslope
    123132        write(num,'(i2.2)') islope
    124133        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    125134    enddo
    126 
    127     write(*,*) "Downloading data for h2o_ice_s_slope done"
     135    write(*,*) "Data for h2o_ice_s downloaded!"
    128136
    129137#ifndef CPP_STD
    130     write(*,*) "Downloading data for watercap_slope ..."
    131138    do islope = 1,nslope
    132139        write(num,'(i2.2)') islope
    133         call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
    134 !        watercap_slope(:,:,:,:) = 0.
    135     enddo
    136     write(*,*) "Downloading data for watercap_slope done"
    137 #endif
    138 
    139     write(*,*) "Downloading data for tsurf_slope ..."
    140 
    141     do islope=1,nslope
     140        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
     141    enddo
     142    write(*,*) "Data for watercap downloaded!"
     143#endif
     144
     145    do islope = 1,nslope
    142146        write(num,'(i2.2)') islope
    143147        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
    144148    enddo
    145 
    146     write(*,*) "Downloading data for tsurf_slope done"
     149    write(*,*) "Data for tsurf downloaded!"
    147150
    148151#ifndef CPP_STD
    149152    if (soil_pem) then
    150         write(*,*) "Downloading data for soiltemp_slope ..."
    151 
    152153        do islope = 1,nslope
    153154            write(num,'(i2.2)') islope
    154155            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
    155156        enddo
    156 
    157         write(*,*) "Downloading data for soiltemp_slope done"
    158         write(*,*) "Downloading data for watersoil_density ..."
     157        write(*,*) "Data for soiltemp downloaded!"
    159158
    160159        do islope = 1,nslope
    161160            write(num,'(i2.2)') islope
    162             call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
    163         enddo
    164 
    165         write(*,*) "Downloading data for  watersoil_density  done"
    166         write(*,*) "Downloading data for  watersurf_density  ..."
     161            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
     162        enddo
     163        write(*,*) "Data for waterdensity_soil downloaded!"
    167164
    168165        do islope = 1,nslope
    169166            write(num,'(i2.2)') islope
    170             call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
    171         enddo
    172 
    173         write(*,*) "Downloading data for  watersurf_density  done"
    174 
     167            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
     168        enddo
     169        write(*,*) "Data for waterdensity_surface downloaded!"
    175170    endif !soil_pem
    176171#endif
    177 else !nslope = 1 no slope, we copy all the values
    178     call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
    179     call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
    180     call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
    181    
    182 #ifndef CPP_STD
    183 !    call get_var3("watercap", watercap_slope(:,:,1,:))
    184     watercap_slope(:,:,1,:) = 0.
    185     watersurf_density_dyn(:,:,:,:) = 0.
    186     watersoil_density_dyn(:,:,:,:,:) = 0.
    187 #endif
    188 
    189     if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
    190 endif !nslope=1
     172endif
    191173
    192174! Compute the minimum over the year for each point
    193 write(*,*) "Computing the min of h2o_ice_slope"
    194 min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
    195 write(*,*) "Computing the min of co2_ice_slope"
     175write(*,*) "Computing the min of h2o_ice..."
     176where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
     177min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap,4)
     178write(*,*) "Computing the min of co2_ice..."
     179where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
    196180min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
    197181
    198 !Compute averages
    199 write(*,*) "Computing average of tsurf"
     182! Compute averages over the year for each point
     183write(*,*) "Computing the average of tsurf..."
    200184tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
    201185
    … …  
    203187    do islope = 1,nslope
    204188        do t = 1,timelen
    205             call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
     189            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
    206190        enddo
    207191    enddo
    … …  
    209193
    210194if (soil_pem) then
    211     write(*,*) "Computing average of tsoil"
     195    write(*,*) "Computing average of tsoil..."
    212196    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    213     write(*,*) "Computing average of watersurf_density"
     197    write(*,*) "Computing average of waterdensity_surface..."
    214198    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
    215199endif
    216200
    217 ! By definition, a density is positive, we get rid of the negative values
    218 where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
    219 where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
    220 
     201! By definition, we get rid of the negative values
     202A = (1./m_co2 - 1./m_noco2)
     203B = 1./m_noco2
    221204do i = 1,iim + 1
    222205    do j = 1,jjm_input + 1
    … …  
    250233            do l = 1,nsoilmx
    251234                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
    252                 do t=1,timelen
     235                do t = 1,timelen
    253236                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
    254237                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
    255238                enddo
    256239            enddo
    257         endif !soil_pem
     240        endif ! soil_pem
    258241        do t = 1,timelen
    259242            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
    … …  
    274257        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
    275258        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
    276     endif !soil_pem
     259    endif ! soil_pem
    277260    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
    278261    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
    … …  
    281264
    282265END SUBROUTINE read_data_PCM
     266
     267!=======================================================================
    283268
    284269SUBROUTINE check_dim(n1,n2,str1,str2)
    … …  
    302287!=======================================================================
    303288
    304 SUBROUTINE get_var1(var,v)
    305 
    306 implicit none
    307 
    308 character(len = *), intent(in)  :: var
     289SUBROUTINE get_var1(filename,v)
     290
     291implicit none
     292
     293character(len = *), intent(in)  :: filename
    309294real, dimension(:), intent(out) :: v
    310295
    311 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
    312 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
     296call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
     297call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
    313298
    314299END SUBROUTINE get_var1
    … …  
    316301!=======================================================================
    317302
    318 SUBROUTINE get_var3(var,v) ! on U grid
    319 
    320 implicit none
    321 
    322 character(len = *),     intent(in)  :: var
     303SUBROUTINE get_var3(filename,v) ! on U grid
     304
     305implicit none
     306
     307character(len = *),     intent(in)  :: filename
    323308real, dimension(:,:,:), intent(out) :: v
    324309
    325 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
    326 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
     310call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
     311call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
    327312
    328313END SUBROUTINE get_var3
    … …  
    330315!=======================================================================
    331316
    332 SUBROUTINE get_var4(var,v)
    333 
    334 implicit none
    335 
    336 character(len = *),       intent(in)  :: var
     317SUBROUTINE get_var4(filename,v)
     318
     319implicit none
     320
     321character(len = *),       intent(in)  :: filename
    337322real, dimension(:,:,:,:), intent(out) :: v
    338323
    339 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
    340 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
     324call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
     325call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
    341326
    342327END SUBROUTINE get_var4

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

    r3097 r3122  
    3232!         Oct 2011 Francois: enable having a 'diagpem.def' file to select
    3333!                            at runtime, which variables to put in file
    34 !         Oct 2023 JB: conversion into Fortran 90 with module for the PEM
     34!         Oct 2023 JBC: conversion into Fortran 90 with module for the PEM
    3535!
    3636!  parametres (input) :
Note: See TracChangeset for help on using the changeset viewer.