Ignore:
Timestamp:
Nov 20, 2024, 3:53:19 PM (25 hours ago)
Author:
jbclement
Message:

PEM:

  • Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
  • Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
  • Cleaning of "glaciers_mod.F90";
  • Optimization for the computation of ice density according to temperature by using a function.

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3345 r3527  
    22
    33implicit none
    4 
    54
    65! Flags for ice management
    … …  
    3332implicit none
    3433
    35 ! arguments
    36 ! ---------
    37 
    3834! Inputs:
    39       integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    40       real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
    41       real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
    42       real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
    43       real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
    44       real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
    45       real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
    46      
     35!--------
     36integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     37real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
     38real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
     39real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
     40real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
     41real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
     42real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
     43
    4744! Ouputs:
    48       real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    49       real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
    50       real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    51 
     45!--------
     46real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     47real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
     48real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    5249! Local
    53       real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
    54       real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    55 
    56 !-----------------------------
     50!------
     51real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
     52real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     53
    5754write(*,*) "Flow of CO2 glaciers"
    58    
    5955call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
    6056call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    … …  
    8379
    8480! Inputs:
    85       integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    86       real,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
    87       real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
     81integer,                       intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     82real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Physical points x Slopes : Distribution of the subgrid slopes
     83real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
     84real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
    8885! Ouputs:
    89       real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    90       real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    91       real,intent(inout) :: flag_h2oflow_mesh(ngrid) ! same but within the mesh
     86real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     87real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
     88real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
    9289! Local
    93       real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    94 
    95 !-----------------------------
     90real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow
     91
    9692write(*,*) "Flow of H2O glaciers"
    97 
    9893call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    9994call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
    … …  
    105100SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    106101
     102use ice_table_mod, only: rho_ice
    107103use abort_pem_mod, only: abort_pem
    108104#ifndef CPP_STD
    … …  
    114110!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    115111!!!
    116 !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
    117 !!!          before initating flow
     112!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
    118113!!!         
    119114!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    … …  
    123118implicit none
    124119
    125 ! arguments
    126 ! --------
    127 
    128120! Inputs
    129       integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
    130       integer,intent(in) :: iflat        ! index of the flat subslope
    131       real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
    132       real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
    133       character(len=3), intent(in) :: name_ice ! Nature of the ice
     121! ------
     122integer,                       intent(in) :: ngrid, nslope  ! # of grid points and subslopes
     123integer,                       intent(in) :: iflat          ! index of the flat subslope
     124real, dimension(nslope),       intent(in) :: def_slope_mean ! Slope field: Values of the subgrid slope angles [deg]
     125real, dimension(ngrid,nslope), intent(in) :: Tice           ! Physical field:  ice temperature [K]
     126character(3),                  intent(in) :: name_ice       ! Nature of ice
    134127! Outputs
    135       real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
     128! -------
     129real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum  thickness before flaw [m]
    136130! Local
    137       DOUBLE PRECISION :: tau_d    ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
    138       real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
    139       integer :: ig,islope ! loop variables
    140       real :: slo_angle
    141 
    142 ! 1. Compute rho
    143     if(name_ice.eq."co2") then
    144       DO ig = 1,ngrid
    145         DO islope = 1,nslope
    146         rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
     131! -----
     132real    :: tau_d      ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
     133integer :: ig, islope ! loop variables
     134real    :: slo_angle
     135
     136select case (trim(adjustl(name_ice)))
     137    case('h2o')
     138        tau_d = 1.e5
     139    case('co2')
    147140        tau_d = 5.e3
    148         ENDDO
    149       ENDDO
    150     elseif (name_ice.eq."h2o") then
    151       DO ig = 1,ngrid
    152         DO islope = 1,nslope
    153           rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
    154           tau_d = 1.e5
    155         ENDDO
    156       ENDDO
    157     else
    158       call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
    159     endif
    160 
    161 ! 3. Compute max thickness
    162     DO ig = 1,ngrid
    163        DO islope = 1,nslope
    164           if(islope.eq.iflat) then
     141    case default
     142        call abort_pem("compute_hmaxglaciers","Type of ice not known!",1)
     143end select
     144
     145do ig = 1,ngrid
     146    do islope = 1,nslope
     147        if (islope == iflat) then
    165148            hmax(ig,islope) = 1.e8
    166           else
     149        else
    167150            slo_angle = abs(def_slope_mean(islope)*pi/180.)
    168             hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
    169           endif
    170        ENDDO
    171     ENDDO
     151            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
     152        endif
     153    enddo
     154enddo
     155
    172156END SUBROUTINE compute_hmaxglaciers
    173157
    … …  
    183167!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    184168
     169use ice_table_mod, only: rho_ice
    185170use abort_pem_mod, only: abort_pem
    186171#ifndef CPP_STD
    … …  
    192177implicit none
    193178
    194 ! arguments
    195 ! --------
    196 
    197179! Inputs
    198       integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
    199       integer, intent(in) :: iflat                      ! index of the flat subslope
    200       real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
    201       real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
    202       real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
    203       real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
    204       character(len=3), intent(in) :: name_ice              ! Nature of the ice
    205 
     180!-------
     181integer,                       intent(in) :: ngrid, nslope  ! # of physical points and subslope
     182integer,                       intent(in) :: iflat          ! index of the flat subslope
     183real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Distribution of the subgrid slopes within the mesh
     184real, dimension(nslope),       intent(in) :: def_slope_mean ! values of the subgrid slopes
     185real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
     186real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
     187character(3),                  intent(in) :: name_ice       ! Nature of the ice
    206188! Outputs
    207       real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
    208       real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
    209       real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
     189!--------
     190real, dimension(ngrid,nslope), intent(inout) :: qice          ! CO2 in the subslope [kg/m^2]
     191real, dimension(ngrid,nslope), intent(inout) :: flag_flow     ! boolean to check if there is flow on a subgrid slope
     192real, dimension(ngrid),        intent(inout) :: flag_flowmesh ! boolean to check if there is flow in the mesh
    210193! Local
    211       integer ig,islope ! loop
    212       real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
    213       integer iaval ! ice will be transfered here
    214 
    215 ! 0. Compute rho
    216     if(name_ice.eq."co2") then
    217       DO ig = 1,ngrid
    218         DO islope = 1,nslope
    219         rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
    220         ENDDO
    221       ENDDO
    222     elseif (name_ice.eq."h2o") then
    223       DO ig = 1,ngrid
    224         DO islope = 1,nslope
    225           rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
    226         ENDDO
    227       ENDDO
    228     else
    229       call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
    230     endif
    231 
    232 ! 1. Compute the transfer of ice
    233 
    234        DO ig = 1,ngrid
    235         DO islope = 1,nslope
    236           IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
     194!------
     195integer :: ig, islope ! loop
     196integer :: iaval      ! ice will be transfered here
     197
     198do ig = 1,ngrid
     199    do islope = 1,nslope
     200        if (islope /= iflat) then ! ice can be infinite on flat ground
    237201! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
    238             IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
    239                   cos(pi*def_slope_mean(islope)/180.)) THEN
     202            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
    240203! Second: determine the flatest slopes possible:
    241                 IF(islope.gt.iflat) THEN
    242                   iaval=islope-1
    243                 ELSE
    244                  iaval=islope+1
    245                 ENDIF
    246                 do while ((iaval.ne.iflat).and.  &
    247                     (subslope_dist(ig,iaval).eq.0))
    248                   IF(iaval.gt.iflat) THEN
    249                      iaval=iaval-1
    250                   ELSE
    251                      iaval=iaval+1
    252                   ENDIF
     204                if (islope > iflat) then
     205                    iaval=islope-1
     206                else
     207                    iaval = islope + 1
     208                endif
     209                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
     210                    if (iaval > iflat) then
     211                        iaval = iaval - 1
     212                    else
     213                        iaval = iaval + 1
     214                    endif
    253215                enddo
    254               qice(ig,iaval) = qice(ig,iaval) + &
    255                (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
    256                cos(pi*def_slope_mean(islope)/180.)) *             &
    257                subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
    258                cos(pi*def_slope_mean(iaval)/180.) /               &
    259                cos(pi*def_slope_mean(islope)/180.)
    260 
    261               qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
    262                cos(pi*def_slope_mean(islope)/180.)
    263 
    264               flag_flow(ig,islope) = 1.
    265               flag_flowmesh(ig) = 1.
    266             ENDIF ! co2ice > hmax
    267           ENDIF ! iflat
    268         ENDDO !islope
    269        ENDDO !ig
     216                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) &
     217                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
     218
     219                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
     220                flag_flow(ig,islope) = 1.
     221                flag_flowmesh(ig) = 1.
     222            endif ! co2ice > hmax
     223        endif ! iflat
     224    enddo !islope
     225enddo !ig
     226
    270227END SUBROUTINE
    271228
    … …  
    281238!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    282239 
    283 use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
     240use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
    284241
    285242implicit none
    … …  
    302259real    :: ave    ! Intermediate to compute average
    303260
    304 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    305261do ig = 1,ngrid
    306262    ave = 0
Note: See TracChangeset for help on using the changeset viewer.