- Timestamp:
- Sep 18, 2024, 4:09:02 PM (3 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3428 r3430 422 422 == 16/09/2024 == JBC 423 423 Corrections and improvements for the launching script. 424 425 == 18/09/2024 == JBC 426 - Correction of the way stopping criteria interact with each other because of a situation where the algorithm does not stop even though a stopping criterion was met. 427 - Update and corrections for the layering algorithm: there is now only one subsurface stratum. 428 - Improvement of the launching script in the case of the 1D PEM: it ends automatically after PCM/PEM run crash. -
trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh
r3428 r3430 175 175 sed -i "s/^k=[0-9]\+$/k=$(echo "3 - $nPCM_ini" | bc -l)/" PCMrun.job 176 176 ./PCMrun.job 177 if [ $? -ne 0 ]; then 178 errlaunch 179 fi 177 180 else # 3D model 178 181 cp PCMrun.job PCMrun${iPCM}.job … … 197 200 sed -i "s/^k=[0-9]\+$/k=$(echo "$i + 2 - $nPCM_ini" | bc -l)/" PCMrun.job 198 201 ./PCMrun.job 202 if [ $? -ne 0 ]; then 203 errlaunch 204 fi 199 205 else # 3D model 200 206 cp PCMrun.job PCMrun${iPCM}.job … … 219 225 if [ $1 -eq 1 ]; then # 1D model 220 226 ./PEMrun.job 227 if [ $? -ne 0 ]; then 228 errlaunch 229 fi 221 230 else # 3D model 222 231 sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job … … 245 254 if [ $1 -eq 1 ]; then # 1D model 246 255 ./PEMrun.job 256 if [ $? -ne 0 ]; then 257 errlaunch 258 fi 247 259 else # 3D model 248 260 sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job … … 351 363 # arg1: model dimension 352 364 relaunchPEM() { 353 iPEM=$ irelaunch365 iPEM=$(($irelaunch + 1)) 354 366 iPCM=$(($nPCM_ini + $nPCM*($irelaunch - 1) + 1)) 355 367 i_myear=$(awk "NR==$(($iPEM + 1)) {print \$1}" "info_PEM.txt") -
trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
r3424 r3430 273 273 ! Creation of three strata at the bottom of the layering 274 274 ! 1) Dry porous deep regolith 275 call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)275 !call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity) 276 276 ! 2) Ice-cemented regolith (ice table) 277 call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)277 !call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.) 278 278 ! 3) Dry porous regolith 279 279 call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity) … … 404 404 405 405 !---- Local variables 406 type(stratum), pointer :: current 407 integer :: ig, islope, k 406 integer :: ig, islope, k 408 407 409 408 !---- Code 410 409 do islope = 1,nslope 411 410 do ig = 1,ngrid 412 current => stratif(ig,islope)%bottom 413 do k = 1,3 414 current%thickness = stratif_array(ig,islope,k,1) 415 current%top_elevation = stratif_array(ig,islope,k,2) 416 current%co2ice_volfrac = stratif_array(ig,islope,k,3) 417 current%h2oice_volfrac = stratif_array(ig,islope,k,4) 418 current%dust_volfrac = stratif_array(ig,islope,k,5) 419 current%air_volfrac = stratif_array(ig,islope,k,6) 420 current => current%up 421 enddo 422 do k = 4,size(stratif_array,3) 411 stratif(ig,islope)%bottom%thickness = stratif_array(ig,islope,1,1) 412 stratif(ig,islope)%bottom%top_elevation = stratif_array(ig,islope,1,2) 413 stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3) 414 stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4) 415 stratif(ig,islope)%bottom%dust_volfrac = stratif_array(ig,islope,1,5) 416 stratif(ig,islope)%bottom%air_volfrac = stratif_array(ig,islope,1,6) 417 do k = 2,size(stratif_array,3) 423 418 if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit 424 419 call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6)) … … 501 496 else ! No, we stop 502 497 write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!' 503 stopPEM = 7498 stopPEM = 8 504 499 exit 505 500 endif … … 507 502 if (h2lift_tot > 0.) then 508 503 write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!' 509 stopPEM = 7504 stopPEM = 8 510 505 endif 511 506 … … 553 548 554 549 ! How much is CO2 ice sublimation inhibited by the top dust lag? 555 currentloc => current1%up556 550 h_toplag = 0. 557 do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 558 h_toplag = h_toplag + currentloc%thickness 559 currentloc => currentloc%up 560 enddo 551 if (associated(current1%up)) then ! If there is an upper stratum 552 currentloc => current1%up 553 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 554 h_toplag = h_toplag + currentloc%thickness 555 if (.not. associated(currentloc%up)) exit 556 currentloc => currentloc%up 557 enddo 558 endif 561 559 h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) 562 560 … … 584 582 if (h2subl_tot > 0.) then 585 583 write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!' 586 stopPEM = 7584 stopPEM = 8 587 585 endif 588 586 ! New stratum at the layering top by condensation of H2O ice … … 609 607 610 608 ! How much is CO2 ice sublimation inhibited by the top dust lag? 611 currentloc => current1%up612 609 h_toplag = 0. 613 do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 614 h_toplag = h_toplag + currentloc%thickness 615 currentloc => currentloc%up 616 enddo 610 if (associated(current1%up)) then ! If there is an upper stratum 611 currentloc => current1%up 612 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 613 h_toplag = h_toplag + currentloc%thickness 614 if (.not. associated(currentloc%up)) exit 615 currentloc => currentloc%up 616 enddo 617 endif 617 618 h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) 618 619 … … 640 641 if (h2subl_tot > 0.) then 641 642 write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!' 642 stopPEM = 7643 stopPEM = 8 643 644 endif 644 645 ! H2O ice sublimation in the considered stratum + New stratum for dust lag … … 649 650 h_dust_old = current2%dust_volfrac*current2%thickness 650 651 651 ! How much is CO2 ice sublimation inhibited by the top dust lag? 652 currentloc => current2%up 652 ! How much is H2O ice sublimation inhibited by the top dust lag? 653 653 h_toplag = 0. 654 do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 655 h_toplag = h_toplag + currentloc%thickness 656 currentloc => currentloc%up 657 enddo 654 if (associated(current2%up)) then ! If there is an upper stratum 655 currentloc => current2%up 656 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol) 657 h_toplag = h_toplag + currentloc%thickness 658 if (.not. associated(currentloc%up)) exit 659 currentloc => currentloc%up 660 enddo 661 endif 658 662 h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) 659 663 … … 681 685 if (h2subl_tot > 0.) then 682 686 write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!' 683 stopPEM = 7687 stopPEM = 8 684 688 endif 685 689 -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3426 r3430 1001 1001 year_iter = year_iter + dt_pem 1002 1002 i_myear = i_myear + dt_pem 1003 if ( year_iter >= year_iter_max) stopPEM = 51004 if ( i_myear >= n_myear) stopPEM = 61003 if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5 1004 if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6 1005 1005 call system_clock(c2) 1006 if ( timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 71007 if (stopPEM >0) then1006 if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7 1007 if (stopPEM <= 0) then 1008 1008 select case (stopPEM) 1009 1009 case(1) … … 1021 1021 case(7) 1022 1022 write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM 1023 case(8) 1024 write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM 1023 1025 case default 1024 1026 write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM -
trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
r3319 r3430 116 116 117 117 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 118 #ifndef CPP_STD 119 ! h2o ice 120 h2o_ice = 0. 121 call get_field("h2o_ice",h2o_ice,found) 122 if (.not. found) then 123 write(*,*)'Pemetat0: failed loading <h2o_ice>' 124 write(*,*)'will reconstruct the values from watercaptag' 125 write(*,*)'with default value ''ini_huge_h2oice''' 126 do ig = 1,ngrid 127 if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice 128 enddo 129 else 130 ! The variations of infinite reservoirs during the PCM years are taken into account 131 h2o_ice = h2o_ice + watercap 132 endif 133 134 ! co2 ice 135 co2_ice = perennial_co2ice 136 #endif 137 138 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 118 139 ! Stratification (layerings) 119 nb_str_max = 3140 nb_str_max = 1 120 141 if (layering_algo) then 121 142 found = inquire_dimension("nb_str_max") … … 149 170 endif 150 171 151 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!152 #ifndef CPP_STD153 ! h2o ice154 h2o_ice = 0.155 call get_field("h2o_ice",h2o_ice,found)156 if (.not. found) then157 write(*,*)'Pemetat0: failed loading <h2o_ice>'158 write(*,*)'will reconstruct the values from watercaptag'159 write(*,*)'with default value ''ini_huge_h2oice'''160 do ig = 1,ngrid161 if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice162 enddo163 else164 ! The variations of infinite reservoirs during the PCM years are taken into account165 h2o_ice = h2o_ice + watercap166 endif167 168 ! co2 ice169 co2_ice = perennial_co2ice170 #endif171 172 172 if (soil_pem) then 173 173 … … 186 186 !!! transition 187 187 delta = depth_breccia 188 TI_PEM(ig,index_breccia +1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &189 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &190 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))191 do iloop=index_breccia +2,index_bedrock188 TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & 189 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + & 190 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) 191 do iloop=index_breccia + 2,index_bedrock 192 192 TI_PEM(ig,iloop,islope) = TI_breccia 193 193 enddo … … 224 224 delta = depth_breccia 225 225 do ig = 1,ngrid 226 if (inertiedat_PEM(ig,index_breccia) .lt.TI_breccia) then226 if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then 227 227 inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & 228 228 (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ & 229 229 ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) 230 230 231 do iloop = index_breccia +2,index_bedrock231 do iloop = index_breccia + 2,index_bedrock 232 232 inertiedat_PEM(ig,iloop) = TI_breccia 233 233 enddo 234 234 else 235 do iloop=index_breccia +1,index_bedrock235 do iloop=index_breccia + 1,index_bedrock 236 236 inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM) 237 237 enddo … … 242 242 delta = depth_bedrock 243 243 do ig = 1,ngrid 244 inertiedat_PEM(ig,index_bedrock +1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &245 (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &246 ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))244 inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ & 245 (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ & 246 ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2)))) 247 247 enddo 248 248 … … 358 358 write(*,*)'There is no "startpem.nc"!' 359 359 360 ! Stratification (layerings)361 nb_str_max = 3362 if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'363 364 360 ! h2o ice 365 361 h2o_ice = 0. … … 372 368 write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.' 373 369 co2_ice = perennial_co2ice 370 371 ! Stratification (layerings) 372 nb_str_max = 1 373 if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.' 374 374 375 375 if (soil_pem) then -
trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90
r3339 r3430 42 42 43 43 !======================================================================= 44 if (stopPEM < 0) return 45 44 46 ! Computation of the present surface of h2o ice still sublimating 45 47 h2oice_now_surf = 0. … … 104 106 105 107 !======================================================================= 108 if (stopPEM < 0) return 109 106 110 ! Computation of the present surface of co2 ice still sublimating 107 111 co2ice_now_surf = 0.
Note: See TracChangeset
for help on using the changeset viewer.