| 1 | SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,& |
|---|
| 2 | ice_depth,TI_PEM) |
|---|
| 3 | #ifndef CPP_STD |
|---|
| 4 | USE comsoil_h, only: inertiedat, volcapa |
|---|
| 5 | USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM |
|---|
| 6 | USE vertical_layers_mod, ONLY: ap,bp |
|---|
| 7 | implicit none |
|---|
| 8 | ! Input: |
|---|
| 9 | INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM |
|---|
| 10 | REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) |
|---|
| 11 | REAL,INTENT(IN) :: p_avg_new |
|---|
| 12 | REAL,INTENT(IN) :: co2ice(ngrid,nslope) |
|---|
| 13 | REAL,INTENT(IN) :: waterice(ngrid,nslope) |
|---|
| 14 | REAL,INTENT(in) :: ice_depth(ngrid,nslope) |
|---|
| 15 | |
|---|
| 16 | ! Outputs: |
|---|
| 17 | |
|---|
| 18 | REAL,INTENT(INOUT) :: TI_PEM(ngrid,nsoil_PEM,nslope) |
|---|
| 19 | |
|---|
| 20 | ! Constants: |
|---|
| 21 | |
|---|
| 22 | REAL :: alpha = 0.2 |
|---|
| 23 | REAL :: beta = 1.08e7 |
|---|
| 24 | REAL :: To = 273.15 |
|---|
| 25 | REAL :: R = 8.314 |
|---|
| 26 | REAL :: L = 51058. |
|---|
| 27 | REAL :: inertie_thresold = 800. ! look for ice |
|---|
| 28 | REAL :: inertie_averaged = 250 ! Mellon et al. 2000 |
|---|
| 29 | REAL :: ice_inertia = 1200 ! Inertia of ice |
|---|
| 30 | REAL :: P610 = 610.0 |
|---|
| 31 | REAL :: m_h2o = 18.01528E-3 |
|---|
| 32 | REAL :: m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) |
|---|
| 33 | REAL :: m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) |
|---|
| 34 | REAL :: A,B,mmean |
|---|
| 35 | |
|---|
| 36 | ! Local variables: |
|---|
| 37 | |
|---|
| 38 | INTEGER :: ig,islope,iloop,iref,k |
|---|
| 39 | REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith |
|---|
| 40 | REAL :: d(ngrid,nsoil_PEM,nslope) |
|---|
| 41 | REAL :: Total_surface |
|---|
| 42 | INTEGER :: ispermanent_co2glaciers(ngrid,nslope) |
|---|
| 43 | INTEGER :: ispermanent_h2oglaciers(ngrid,nslope) |
|---|
| 44 | |
|---|
| 45 | ! 0. Initialisation |
|---|
| 46 | |
|---|
| 47 | ! do ig = 1,ngrid |
|---|
| 48 | ! do islope = 1,nslope |
|---|
| 49 | ! if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then |
|---|
| 50 | ! ispermanent_h2oglaciers(ig,islope) = 1 |
|---|
| 51 | ! else |
|---|
| 52 | ! ispermanent_h2oglaciers(ig,islope) = 0 |
|---|
| 53 | ! endif |
|---|
| 54 | ! |
|---|
| 55 | ! if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then |
|---|
| 56 | ! ispermanent_co2glaciers(ig,islope) = 1 |
|---|
| 57 | ! else |
|---|
| 58 | ! ispermanent_co2glaciers(ig,islope) = 0 |
|---|
| 59 | ! endif |
|---|
| 60 | ! enddo |
|---|
| 61 | ! enddo |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | ! 1.Ice TI feedback |
|---|
| 66 | |
|---|
| 67 | ! do islope = 1,nslope |
|---|
| 68 | ! call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope), TI_PEM(:,:,islope)) |
|---|
| 69 | ! enddo |
|---|
| 70 | |
|---|
| 71 | ! 2. Modification of the regolith thermal inertia. |
|---|
| 72 | |
|---|
| 73 | do ig=1,ngrid |
|---|
| 74 | do islope=1,nslope |
|---|
| 75 | do iloop =1,n_1km |
|---|
| 76 | d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta))) |
|---|
| 77 | if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then ! we are modifying the regolith properties, not ice |
|---|
| 78 | TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta))) |
|---|
| 79 | |
|---|
| 80 | endif |
|---|
| 81 | enddo |
|---|
| 82 | enddo |
|---|
| 83 | enddo |
|---|
| 84 | |
|---|
| 85 | ! 3. Build new TI for the PEM |
|---|
| 86 | ! a) For the regolith |
|---|
| 87 | do ig=1,ngrid |
|---|
| 88 | do islope=1,nslope |
|---|
| 89 | |
|---|
| 90 | if (ice_depth(ig,islope).gt. -1.e-10) then |
|---|
| 91 | ! 3.0 FIrst if permanent ice |
|---|
| 92 | if (ice_depth(ig,islope).lt. 1e-10) then |
|---|
| 93 | do iloop = 1,nsoil_PEM |
|---|
| 94 | TI_PEM(ig,iloop,islope)=max(ice_inertia,inertiedat_PEM(ig,iloop)) |
|---|
| 95 | enddo |
|---|
| 96 | else |
|---|
| 97 | ! 4.1 find the index of the mixed layer |
|---|
| 98 | iref=0 ! initialize iref |
|---|
| 99 | do k=1,nsoil_PEM ! loop on layers |
|---|
| 100 | if (ice_depth(ig,islope).ge.layer_PEM(k)) then |
|---|
| 101 | iref=k ! pure regolith layer up to here |
|---|
| 102 | else |
|---|
| 103 | ! correct iref was obtained in previous cycle |
|---|
| 104 | exit |
|---|
| 105 | endif |
|---|
| 106 | enddo |
|---|
| 107 | |
|---|
| 108 | ! 4.2 Build the new ti |
|---|
| 109 | do iloop=1,iref |
|---|
| 110 | TI_PEM(ig,iloop,islope) =TI_PEM(ig,iloop,islope) |
|---|
| 111 | enddo |
|---|
| 112 | if (iref.lt.nsoil_PEM) then |
|---|
| 113 | if (iref.ne.0) then |
|---|
| 114 | ! mixed layer |
|---|
| 115 | TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ & |
|---|
| 116 | (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ & |
|---|
| 117 | ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2)))) |
|---|
| 118 | else ! first layer is already a mixed layer |
|---|
| 119 | ! (ie: take layer(iref=0)=0) |
|---|
| 120 | TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ & |
|---|
| 121 | (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ & |
|---|
| 122 | ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2)))) |
|---|
| 123 | endif ! of if (iref.ne.0) |
|---|
| 124 | ! lower layers of pure ice |
|---|
| 125 | do iloop=iref+2,nsoil_PEM |
|---|
| 126 | TI_PEM(ig,iloop,islope)=ice_inertia |
|---|
| 127 | enddo |
|---|
| 128 | endif ! of if (iref.lt.(nsoilmx)) |
|---|
| 129 | endif ! permanent glaciers |
|---|
| 130 | endif ! depth > 0 |
|---|
| 131 | enddo !islope |
|---|
| 132 | enddo !ig |
|---|
| 133 | |
|---|
| 134 | !======================================================================= |
|---|
| 135 | RETURN |
|---|
| 136 | #endif |
|---|
| 137 | END |
|---|
