Changeset 699 for trunk/LMDZ.GENERIC/libf
- Timestamp:
- Jun 6, 2012, 12:47:56 PM (13 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf
- Files:
-
- 3 edited
-
dyn3d/calfis.F (modified) (1 diff)
-
dyn3d/newstart.F (modified) (27 diffs)
-
phystd/phyetat0.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/dyn3d/calfis.F
r135 r699 178 178 latfi(ngridmx)= rlatu(jjp1) 179 179 lonfi(ngridmx)= 0. 180 181 ! build airefi(), mesh area on physics grid 180 182 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) 183 ! Poles are single points on physics grid 184 airefi(1)=airefi(1)*iim 185 airefi(ngridmx)=airefi(ngridmx)*iim 181 186 182 187 CALL inifis(ngridmx,llm,day_ini,daysec,dtphys, -
trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F
r649 r699 145 145 character(len=20) :: txt ! to store some text 146 146 integer :: count 147 148 ! MONS data: 149 real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O 150 real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2) 151 ! coefficient to apply to convert d21 to 'true' depth (m) 152 real :: MONS_coeff 153 real :: MONS_coeffS ! coeff for southern hemisphere 154 real :: MONS_coeffN ! coeff for northern hemisphere 155 ! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth 147 real :: profile(llm+1) ! to store an atmospheric profile + surface value 156 148 157 149 ! added by RW for test … … 404 396 latfi(ngridmx)=rlatu(jjp1) 405 397 lonfi(ngridmx)=0. 398 399 ! build airefi(), mesh area on physics grid 406 400 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) 401 ! Poles are single points on physics grid 402 airefi(1)=airefi(1)*iim 403 airefi(ngridmx)=airefi(ngridmx)*iim 407 404 408 405 c======================================================================= … … 479 476 write(*,*) 'q=0 : ALL tracer =zero' 480 477 write(*,*) 'q=x : give a specific uniform value to one tracer' 481 write(*,*) 'ini_q : tracers initialisation for chemistry, water an 482 $d ice ' 483 write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and 484 $ice ' 485 write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on 486 $ly ' 478 write(*,*) 'q=profile : specify a profile for a tracer' 479 ! write(*,*) 'ini_q : tracers initialisation for chemistry, water an 480 ! $d ice ' 481 ! write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and 482 ! $ice ' 483 ! write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on 484 ! $ly ' 487 485 write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' 488 486 write(*,*) 'watercapn : H20 ice on permanent N polar cap ' … … 502 500 write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur 503 501 &face values' 504 ! write(*,*) 'subsoilice_n : Put deep underground ice layer in northe505 ! &rn hemisphere'506 ! write(*,*) 'subsoilice_s : Put deep underground ice layer in southe507 ! &rn hemisphere'508 ! write(*,*) 'mons_ice : Put underground ice layer according to MONS-509 ! &derived data'510 502 511 503 write(*,*) … … 522 514 c 'flat : no topography ("aquaplanet")' 523 515 c ------------------------------------- 524 if ( modif(1:len_trim(modif)) .eq. 'flat') then516 if (trim(modif) .eq. 'flat') then 525 517 c set topo to zero 526 CALL initial0(ip1jmp1,z_reel)527 CALL multscal(ip1jmp1,z_reel,g,phis)518 z_reel(1:iip1,1:jjp1)=0 519 phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1) 528 520 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) 529 521 write(*,*) 'topography set to zero.' … … 556 548 c 'nuketharsis : no tharsis bulge for Early Mars' 557 549 c --------------------------------------------- 558 else if ( modif(1:len_trim(modif)) .eq. 'nuketharsis') then550 else if (trim(modif) .eq. 'nuketharsis') then 559 551 560 552 DO j=1,jjp1 … … 583 575 c bilball : uniform albedo, thermal inertia 584 576 c ----------------------------------------- 585 else if ( modif(1:len_trim(modif)) .eq. 'bilball') then577 else if (trim(modif) .eq. 'bilball') then 586 578 write(*,*) 'constante albedo and iner.therm:' 587 579 write(*,*) 'New value for albedo (ex: 0.25) ?' … … 610 602 c coldspole : sous-sol de la calotte sud toujours froid 611 603 c ----------------------------------------------------- 612 else if ( modif(1:len_trim(modif)) .eq. 'coldspole') then604 else if (trim(modif) .eq. 'coldspole') then 613 605 write(*,*)'new value for the subsurface temperature', 614 606 & ' beneath the permanent southern polar cap ? (eg: 141 K)' … … 661 653 c ptot : Modification of the total pressure: ice + current atmosphere 662 654 c ------------------------------------------------------------------- 663 else if ( modif(1:len_trim(modif)).eq.'ptot') then655 else if (trim(modif).eq.'ptot') then 664 656 665 657 ! check if we have a co2_ice surface tracer: … … 789 781 c q=0 : set tracers to zero 790 782 c ------------------------- 791 else if ( modif(1:len_trim(modif)).eq.'q=0') then783 else if (trim(modif).eq.'q=0') then 792 784 c mise a 0 des q (traceurs) 793 785 write(*,*) 'Tracers set to 0 (1.E-30 in fact)' … … 811 803 c q=x : initialise tracer manually 812 804 c -------------------------------- 813 else if ( modif(1:len_trim(modif)).eq.'q=x') then805 else if (trim(modif).eq.'q=x') then 814 806 write(*,*) 'Which tracer do you want to modify ?' 815 807 do iq=1,nqmx … … 835 827 ENDDO 836 828 837 c ini_q : Initialize tracers for chemistry 838 c ----------------------------------------------- 839 else if (modif(1:len_trim(modif)) .eq. 'ini_q') then 840 c For more than 32 layers, possible to initiate thermosphere only 841 thermo=0 842 yes=' ' 843 if (llm.gt.32) then 844 do while ((yes.ne.'y').and.(yes.ne.'n')) 845 write(*,*)'', 846 & 'initialisation for thermosphere only? (y/n)' 847 read(*,fmt='(a)') yes 848 if (yes.eq.'y') then 849 thermo=1 850 else 851 thermo=0 852 endif 853 enddo 854 endif 855 856 c call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi, 857 c $ thermo,qsurf) 858 write(*,*) 'Chemical species initialized' 859 860 if (thermo.eq.0) then 861 c mise a 0 des qsurf (traceurs a la surface) 862 DO iq =1, nqmx 863 DO ig=1,ngridmx 864 qsurf(ig,iq)=0. 865 ENDDO 866 ENDDO 867 endif 868 869 c ini_q-H2O : as above except for the water vapour tracer 870 c ------------------------------------------------------ 871 else if (modif(1:len_trim(modif)) .eq. 'ini_q-H2O') then 872 ! for more than 32 layers, possible to initiate thermosphere only 873 thermo=0 874 yes=' ' 875 if(llm.gt.32) then 876 do while ((yes.ne.'y').and.(yes.ne.'n')) 877 write(*,*)'', 878 & 'initialisation for thermosphere only? (y/n)' 879 read(*,fmt='(a)') yes 880 if (yes.eq.'y') then 881 thermo=1 882 else 883 thermo=0 884 endif 885 enddo 886 endif 887 c call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi, 888 c $ thermo,qsurf) 889 c write(*,*) 'Initialized chem. species exept last (H2O)' 890 891 if (thermo.eq.0) then 892 c set surface tracers to zero, except water ice 893 DO iq =1, nqmx 894 if (iq.ne.igcm_h2o_ice) then 895 DO ig=1,ngridmx 896 qsurf(ig,iq)=0. 897 ENDDO 898 endif 899 ENDDO 900 endif 901 902 c ini_q-iceH2O : as above exept for ice et H2O 903 c ----------------------------------------------- 904 else if (modif(1:len_trim(modif)) .eq. 'ini_q-iceH2O') then 905 c For more than 32 layers, possible to initiate thermosphere only 906 thermo=0 907 yes=' ' 908 if(llm.gt.32) then 909 do while ((yes.ne.'y').and.(yes.ne.'n')) 910 write(*,*)'', 911 & 'initialisation for thermosphere only? (y/n)' 912 read(*,fmt='(a)') yes 913 if (yes.eq.'y') then 914 thermo=1 915 else 916 thermo=0 917 endif 918 enddo 919 endif 920 921 c call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi, 922 c $ thermo,qsurf) 923 c write(*,*) 'Initialized chem. species exept ice and H2O' 924 925 if (thermo.eq.0) then 926 c set surface tracers to zero, except water ice 927 DO iq =1, nqmx 928 if (iq.ne.igcm_h2o_ice) then 929 DO ig=1,ngridmx 930 qsurf(ig,iq)=0. 931 ENDDO 932 endif 933 ENDDO 934 endif 829 c q=profile : initialize tracer with a given profile 830 c -------------------------------------------------- 831 else if (trim(modif) .eq. 'q=profile') then 832 write(*,*) 'Tracer profile will be sought in ASCII file' 833 write(*,*) "'profile_tracer' where 'tracer' is tracer name" 834 write(*,*) "(one value per line in file; starting with" 835 write(*,*) "surface value, the 1st atmospheric layer" 836 write(*,*) "followed by 2nd, etc. up to top of atmosphere)" 837 write(*,*) 'Which tracer do you want to set?' 838 do iq=1,nqmx 839 write(*,*)iq,' : ',trim(tnom(iq)) 840 enddo 841 write(*,*) '(choose between 1 and ',nqmx,')' 842 read(*,*) iq 843 if ((iq.lt.1).or.(iq.gt.nqmx)) then 844 ! wrong value for iq, go back to menu 845 write(*,*) "wrong input value:",iq 846 cycle 847 endif 848 ! look for input file 'profile_tracer' 849 txt="profile_"//trim(tnom(iq)) 850 open(41,file=trim(txt),status='old',form='formatted', 851 & iostat=ierr) 852 if (ierr.eq.0) then 853 ! OK, found file 'profile_...', load the profile 854 do l=1,llm+1 855 read(41,*,iostat=ierr) profile(l) 856 if (ierr.ne.0) then ! something went wrong 857 exit ! quit loop 858 endif 859 enddo 860 if (ierr.eq.0) then 861 ! initialize tracer values 862 qsurf(:,iq)=profile(1) 863 do l=1,llm 864 q(:,:,l,iq)=profile(l+1) 865 enddo 866 write(*,*)'OK, tracer ',trim(tnom(iq)),' initialized ', 867 & 'using values from file ',trim(txt) 868 else 869 write(*,*)'problem reading file ',trim(txt),' !' 870 write(*,*)'No modifications to tracer ',trim(tnom(iq)) 871 endif 872 else 873 write(*,*)'Could not find file ',trim(txt),' !' 874 write(*,*)'No modifications to tracer ',trim(tnom(iq)) 875 endif 876 935 877 936 878 c wetstart : wet atmosphere with a north to south gradient 937 879 c -------------------------------------------------------- 938 else if ( modif(1:len_trim(modif)) .eq. 'wetstart') then880 else if (trim(modif) .eq. 'wetstart') then 939 881 ! check that there is indeed a water vapor tracer 940 882 if (igcm_h2o_vap.eq.0) then … … 957 899 c noglacier : remove tropical water ice (to initialize high res sim) 958 900 c -------------------------------------------------- 959 else if ( modif(1:len_trim(modif)) .eq. 'noglacier') then901 else if (trim(modif) .eq. 'noglacier') then 960 902 if (igcm_h2o_ice.eq.0) then 961 903 write(*,*) "No water ice tracer! Can't use this option" … … 974 916 c watercapn : H20 ice on permanent northern cap 975 917 c -------------------------------------------------- 976 else if ( modif(1:len_trim(modif)) .eq. 'watercapn') then918 else if (trim(modif) .eq. 'watercapn') then 977 919 if (igcm_h2o_ice.eq.0) then 978 920 write(*,*) "No water ice tracer! Can't use this option" … … 1016 958 c watercaps : H20 ice on permanent southern cap 1017 959 c ------------------------------------------------- 1018 else if ( modif(1:len_trim(modif)) .eq. 'watercaps') then960 else if (trim(modif) .eq. 'watercaps') then 1019 961 if (igcm_h2o_ice.eq.0) then 1020 962 write(*,*) "No water ice tracer! Can't use this option" … … 1057 999 c noacglac : H2O ice across highest terrain 1058 1000 c -------------------------------------------- 1059 else if ( modif(1:len_trim(modif)) .eq. 'noacglac') then1001 else if (trim(modif) .eq. 'noacglac') then 1060 1002 if (igcm_h2o_ice.eq.0) then 1061 1003 write(*,*) "No water ice tracer! Can't use this option" … … 1086 1028 c oborealis : H2O oceans across Vastitas Borealis 1087 1029 c ----------------------------------------------- 1088 else if ( modif(1:len_trim(modif)) .eq. 'oborealis') then1030 else if (trim(modif) .eq. 'oborealis') then 1089 1031 if (igcm_h2o_ice.eq.0) then 1090 1032 write(*,*) "No water ice tracer! Can't use this option" … … 1118 1060 c iborealis : H2O ice in Northern plains 1119 1061 c -------------------------------------- 1120 else if ( modif(1:len_trim(modif)) .eq. 'iborealis') then1062 else if (trim(modif) .eq. 'iborealis') then 1121 1063 if (igcm_h2o_ice.eq.0) then 1122 1064 write(*,*) "No water ice tracer! Can't use this option" … … 1142 1084 c oceanball : H2O liquid everywhere 1143 1085 c ---------------------------- 1144 else if ( modif(1:len_trim(modif)) .eq. 'oceanball') then1086 else if (trim(modif) .eq. 'oceanball') then 1145 1087 if (igcm_h2o_ice.eq.0) then 1146 1088 write(*,*) "No water ice tracer! Can't use this option" … … 1171 1113 c iceball : H2O ice everywhere 1172 1114 c ---------------------------- 1173 else if ( modif(1:len_trim(modif)) .eq. 'iceball') then1115 else if (trim(modif) .eq. 'iceball') then 1174 1116 if (igcm_h2o_ice.eq.0) then 1175 1117 write(*,*) "No water ice tracer! Can't use this option" … … 1198 1140 c isotherm : Isothermal temperatures and no winds 1199 1141 c ----------------------------------------------- 1200 else if ( modif(1:len_trim(modif)) .eq. 'isotherm') then1142 else if (trim(modif) .eq. 'isotherm') then 1201 1143 1202 1144 write(*,*)'Isothermal temperature of the atmosphere, … … 1206 1148 if(ierr.ne.0) goto 203 1207 1149 1208 do ig=1, ngridmx 1209 tsurf(ig) = Tiso 1210 end do 1211 do l=2,nsoilmx 1212 do ig=1, ngridmx 1213 tsoil(ig,l) = Tiso 1214 end do 1215 end do 1216 DO j=1,jjp1 1217 DO i=1,iim 1218 Do l=1,llm 1219 Tset(i,j,l)=Tiso 1220 end do 1221 end do 1222 end do 1150 tsurf(1:ngridmx)=Tiso 1151 1152 tsoil(1:ngridmx,1:nsoilmx)=Tiso 1153 1154 Tset(1:iip1,1:jjp1,1:llm)=Tiso 1223 1155 flagtset=.true. 1224 call initial0(llm*ip1jmp1,ucov) 1225 call initial0(llm*ip1jm,vcov) 1226 call initial0(ngridmx*(llm+1),q2) 1156 1157 ucov(1:iip1,1:jjp1,1:llm)=0 1158 vcov(1:iip1,1:jjm,1:llm)=0 1159 q2(1:ngridmx,1:nlayermx+1)=0 1227 1160 1228 1161 c radequi : Radiative equilibrium profile of temperatures and no winds 1229 1162 c -------------------------------------------------------------------- 1230 else if ( modif(1:len_trim(modif)) .eq. 'radequi') then1163 else if (trim(modif) .eq. 'radequi') then 1231 1164 1232 1165 write(*,*)'radiative equilibrium temperature profile' … … 1252 1185 end do 1253 1186 flagtset=.true. 1254 call initial0(llm*ip1jmp1,ucov)1255 call initial0(llm*ip1jm,vcov)1256 call initial0(ngridmx*(llm+1),q2)1187 ucov(1:iip1,1:jjp1,1:llm)=0 1188 vcov(1:iip1,1:jjm,1:llm)=0 1189 q2(1:ngridmx,1:nlayermx+1)=0 1257 1190 1258 1191 c coldstart : T set 1K above CO2 frost point and no winds 1259 1192 c ------------------------------------------------ 1260 else if ( modif(1:len_trim(modif)) .eq. 'coldstart') then1193 else if (trim(modif) .eq. 'coldstart') then 1261 1194 1262 1195 write(*,*)'set temperature of the atmosphere,' … … 1288 1221 1289 1222 flagtset=.true. 1290 call initial0(llm*ip1jmp1,ucov)1291 call initial0(llm*ip1jm,vcov)1292 call initial0(ngridmx*(llm+1),q2)1223 ucov(1:iip1,1:jjp1,1:llm)=0 1224 vcov(1:iip1,1:jjm,1:llm)=0 1225 q2(1:ngridmx,1:nlayermx+1)=0 1293 1226 1294 1227 1295 1228 c co2ice=0 : remove CO2 polar ice caps' 1296 1229 c ------------------------------------------------ 1297 else if ( modif(1:len_trim(modif)) .eq. 'co2ice=0') then1230 else if (trim(modif) .eq. 'co2ice=0') then 1298 1231 ! check that there is indeed a co2_ice tracer ... 1299 1232 if (igcm_co2_ice.ne.0) then … … 1310 1243 ! ---------------------------------------------------------------------- 1311 1244 1312 else if ( modif(1:len_trim(modif)).eq.'therm_ini_s') then1245 else if (trim(modif).eq.'therm_ini_s') then 1313 1246 ! write(*,*)"surfithfi(1):",surfithfi(1) 1314 1247 do isoil=1,nsoilmx … … 1331 1264 1332 1265 1333 1334 1335 c$$$! subsoilice_n: Put deep ice layer in northern hemisphere soil 1336 c$$$! ------------------------------------------------------------ 1337 c$$$ 1338 c$$$ else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then 1339 c$$$ 1340 c$$$ write(*,*)'From which latitude (in deg.), up to the north pole, 1341 c$$$ &should we put subterranean ice?' 1342 c$$$ ierr=1 1343 c$$$ do while (ierr.ne.0) 1344 c$$$ read(*,*,iostat=ierr) val 1345 c$$$ if (ierr.eq.0) then ! got a value 1346 c$$$ ! do a sanity check 1347 c$$$ if((val.lt.0.).or.(val.gt.90)) then 1348 c$$$ write(*,*)'Latitude should be between 0 and 90 deg. !!!' 1349 c$$$ ierr=1 1350 c$$$ else ! find corresponding jref (nearest latitude) 1351 c$$$ ! note: rlatu(:) contains decreasing values of latitude 1352 c$$$ ! starting from PI/2 to -PI/2 1353 c$$$ do j=1,jjp1 1354 c$$$ if ((rlatu(j)*180./pi.ge.val).and. 1355 c$$$ & (rlatu(j+1)*180./pi.le.val)) then 1356 c$$$ ! find which grid point is nearest to val: 1357 c$$$ if (abs(rlatu(j)*180./pi-val).le. 1358 c$$$ & abs((rlatu(j+1)*180./pi-val))) then 1359 c$$$ jref=j 1360 c$$$ else 1361 c$$$ jref=j+1 1362 c$$$ endif 1363 c$$$ 1364 c$$$ write(*,*)'Will use nearest grid latitude which is:', 1365 c$$$ & rlatu(jref)*180./pi 1366 c$$$ endif 1367 c$$$ enddo ! of do j=1,jjp1 1368 c$$$ endif ! of if((val.lt.0.).or.(val.gt.90)) 1369 c$$$ endif !of if (ierr.eq.0) 1370 c$$$ enddo ! of do while 1371 c$$$ 1372 c$$$ ! Build layers() (as in soil_settings.F) 1373 c$$$ val2=sqrt(mlayer(0)*mlayer(1)) 1374 c$$$ val3=mlayer(1)/mlayer(0) 1375 c$$$ do isoil=1,nsoilmx 1376 c$$$ layer(isoil)=val2*(val3**(isoil-1)) 1377 c$$$ enddo 1378 c$$$ 1379 c$$$ write(*,*)'At which depth (in m.) does the ice layer begin?' 1380 c$$$ write(*,*)'(currently, the deepest soil layer extends down to:' 1381 c$$$ & ,layer(nsoilmx),')' 1382 c$$$ ierr=1 1383 c$$$ do while (ierr.ne.0) 1384 c$$$ read(*,*,iostat=ierr) val2 1385 c$$$! write(*,*)'val2:',val2,'ierr=',ierr 1386 c$$$ if (ierr.eq.0) then ! got a value, but do a sanity check 1387 c$$$ if(val2.gt.layer(nsoilmx)) then 1388 c$$$ write(*,*)'Depth should be less than ',layer(nsoilmx) 1389 c$$$ ierr=1 1390 c$$$ endif 1391 c$$$ if(val2.lt.layer(1)) then 1392 c$$$ write(*,*)'Depth should be more than ',layer(1) 1393 c$$$ ierr=1 1394 c$$$ endif 1395 c$$$ endif 1396 c$$$ enddo ! of do while 1397 c$$$ 1398 c$$$ ! find the reference index iref the depth corresponds to 1399 c$$$! if (val2.lt.layer(1)) then 1400 c$$$! iref=1 1401 c$$$! else 1402 c$$$ do isoil=1,nsoilmx-1 1403 c$$$ if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 1404 c$$$ & then 1405 c$$$ iref=isoil 1406 c$$$ exit 1407 c$$$ endif 1408 c$$$ enddo 1409 c$$$! endif 1410 c$$$ 1411 c$$$! write(*,*)'iref:',iref,' jref:',jref 1412 c$$$! write(*,*)'layer',layer 1413 c$$$! write(*,*)'mlayer',mlayer 1414 c$$$ 1415 c$$$ ! thermal inertia of the ice: 1416 c$$$ ierr=1 1417 c$$$ do while (ierr.ne.0) 1418 c$$$ write(*,*)'What is the value of subterranean ice thermal inert 1419 c$$$ &ia? (e.g.: 2000)' 1420 c$$$ read(*,*,iostat=ierr)iceith 1421 c$$$ enddo ! of do while 1422 c$$$ 1423 c$$$ ! recast ithfi() onto ith() 1424 c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1425 c$$$ 1426 c$$$ do j=1,jref 1427 c$$$! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi 1428 c$$$ do i=1,iip1 ! loop on longitudes 1429 c$$$ ! Build "equivalent" thermal inertia for the mixed layer 1430 c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ 1431 c$$$ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ 1432 c$$$ & ((layer(iref+1)-val2)/(iceith)**2))) 1433 c$$$ ! Set thermal inertia of lower layers 1434 c$$$ do isoil=iref+2,nsoilmx 1435 c$$$ ith(i,j,isoil)=iceith ! ice 1436 c$$$ enddo 1437 c$$$ enddo ! of do i=1,iip1 1438 c$$$ enddo ! of do j=1,jjp1 1439 c$$$ 1440 c$$$ 1441 c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) 1442 c$$$ 1443 c$$$! do i=1,nsoilmx 1444 c$$$! write(*,*)'i:',i,'ithfi(1,i):',ithfi(1,i) 1445 c$$$! enddo 1446 1447 1448 c$$$! subsoilice_s: Put deep ice layer in southern hemisphere soil 1449 c$$$! ------------------------------------------------------------ 1450 c$$$ 1451 c$$$ else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then 1452 c$$$ 1453 c$$$ write(*,*)'From which latitude (in deg.), down to the south pol 1454 c$$$ &e, should we put subterranean ice?' 1455 c$$$ ierr=1 1456 c$$$ do while (ierr.ne.0) 1457 c$$$ read(*,*,iostat=ierr) val 1458 c$$$ if (ierr.eq.0) then ! got a value 1459 c$$$ ! do a sanity check 1460 c$$$ if((val.gt.0.).or.(val.lt.-90)) then 1461 c$$$ write(*,*)'Latitude should be between 0 and -90 deg. !!!' 1462 c$$$ ierr=1 1463 c$$$ else ! find corresponding jref (nearest latitude) 1464 c$$$ ! note: rlatu(:) contains decreasing values of latitude 1465 c$$$ ! starting from PI/2 to -PI/2 1466 c$$$ do j=1,jjp1 1467 c$$$ if ((rlatu(j)*180./pi.ge.val).and. 1468 c$$$ & (rlatu(j+1)*180./pi.le.val)) then 1469 c$$$ ! find which grid point is nearest to val: 1470 c$$$ if (abs(rlatu(j)*180./pi-val).le. 1471 c$$$ & abs((rlatu(j+1)*180./pi-val))) then 1472 c$$$ jref=j 1473 c$$$ else 1474 c$$$ jref=j+1 1475 c$$$ endif 1476 c$$$ 1477 c$$$ write(*,*)'Will use nearest grid latitude which is:', 1478 c$$$ & rlatu(jref)*180./pi 1479 c$$$ endif 1480 c$$$ enddo ! of do j=1,jjp1 1481 c$$$ endif ! of if((val.lt.0.).or.(val.gt.90)) 1482 c$$$ endif !of if (ierr.eq.0) 1483 c$$$ enddo ! of do while 1484 c$$$ 1485 c$$$ ! Build layers() (as in soil_settings.F) 1486 c$$$ val2=sqrt(mlayer(0)*mlayer(1)) 1487 c$$$ val3=mlayer(1)/mlayer(0) 1488 c$$$ do isoil=1,nsoilmx 1489 c$$$ layer(isoil)=val2*(val3**(isoil-1)) 1490 c$$$ enddo 1491 c$$$ 1492 c$$$ write(*,*)'At which depth (in m.) does the ice layer begin?' 1493 c$$$ write(*,*)'(currently, the deepest soil layer extends down to:' 1494 c$$$ & ,layer(nsoilmx),')' 1495 c$$$ ierr=1 1496 c$$$ do while (ierr.ne.0) 1497 c$$$ read(*,*,iostat=ierr) val2 1498 c$$$! write(*,*)'val2:',val2,'ierr=',ierr 1499 c$$$ if (ierr.eq.0) then ! got a value, but do a sanity check 1500 c$$$ if(val2.gt.layer(nsoilmx)) then 1501 c$$$ write(*,*)'Depth should be less than ',layer(nsoilmx) 1502 c$$$ ierr=1 1503 c$$$ endif 1504 c$$$ if(val2.lt.layer(1)) then 1505 c$$$ write(*,*)'Depth should be more than ',layer(1) 1506 c$$$ ierr=1 1507 c$$$ endif 1508 c$$$ endif 1509 c$$$ enddo ! of do while 1510 c$$$ 1511 c$$$ ! find the reference index iref the depth corresponds to 1512 c$$$ do isoil=1,nsoilmx-1 1513 c$$$ if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 1514 c$$$ & then 1515 c$$$ iref=isoil 1516 c$$$ exit 1517 c$$$ endif 1518 c$$$ enddo 1519 c$$$ 1520 c$$$! write(*,*)'iref:',iref,' jref:',jref 1521 c$$$ 1522 c$$$ ! thermal inertia of the ice: 1523 c$$$ ierr=1 1524 c$$$ do while (ierr.ne.0) 1525 c$$$ write(*,*)'What is the value of subterranean ice thermal inert 1526 c$$$ &ia? (e.g.: 2000)' 1527 c$$$ read(*,*,iostat=ierr)iceith 1528 c$$$ enddo ! of do while 1529 c$$$ 1530 c$$$ ! recast ithfi() onto ith() 1531 c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1532 c$$$ 1533 c$$$ do j=jref,jjp1 1534 c$$$! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi 1535 c$$$ do i=1,iip1 ! loop on longitudes 1536 c$$$ ! Build "equivalent" thermal inertia for the mixed layer 1537 c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ 1538 c$$$ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ 1539 c$$$ & ((layer(iref+1)-val2)/(iceith)**2))) 1540 c$$$ ! Set thermal inertia of lower layers 1541 c$$$ do isoil=iref+2,nsoilmx 1542 c$$$ ith(i,j,isoil)=iceith ! ice 1543 c$$$ enddo 1544 c$$$ enddo ! of do i=1,iip1 1545 c$$$ enddo ! of do j=jref,jjp1 1546 c$$$ 1547 c$$$ 1548 c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) 1549 1550 1551 c$$$c 'mons_ice' : use MONS data to build subsurface ice table 1552 c$$$c -------------------------------------------------------- 1553 c$$$ else if (modif(1:len_trim(modif)).eq.'mons_ice') then 1554 c$$$ 1555 c$$$ ! 1. Load MONS data 1556 c$$$ call load_MONS_data(MONS_Hdn,MONS_d21) 1557 c$$$ 1558 c$$$ ! 2. Get parameters from user 1559 c$$$ ierr=1 1560 c$$$ do while (ierr.ne.0) 1561 c$$$ write(*,*) "Coefficient to apply to MONS 'depth' in Northern", 1562 c$$$ & " Hemisphere?" 1563 c$$$ write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)" 1564 c$$$ read(*,*,iostat=ierr) MONS_coeffN 1565 c$$$ enddo 1566 c$$$ ierr=1 1567 c$$$ do while (ierr.ne.0) 1568 c$$$ write(*,*) "Coefficient to apply to MONS 'depth' in Southern", 1569 c$$$ & " Hemisphere?" 1570 c$$$ write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)" 1571 c$$$ read(*,*,iostat=ierr) MONS_coeffS 1572 c$$$ enddo 1573 c$$$ ierr=1 1574 c$$$ do while (ierr.ne.0) 1575 c$$$ write(*,*) "Value of subterranean ice thermal inertia?" 1576 c$$$ write(*,*) " (e.g.: 2000, or perhaps 2290)" 1577 c$$$ read(*,*,iostat=ierr) iceith 1578 c$$$ enddo 1579 c$$$ 1580 c$$$ ! 3. Build subterranean thermal inertia 1581 c$$$ 1582 c$$$ ! initialise subsurface inertia with reference surface values 1583 c$$$ do isoil=1,nsoilmx 1584 c$$$ ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx) 1585 c$$$ enddo 1586 c$$$ ! recast ithfi() onto ith() 1587 c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1588 c$$$ 1589 c$$$ do i=1,iip1 ! loop on longitudes 1590 c$$$ do j=1,jjp1 ! loop on latitudes 1591 c$$$ ! set MONS_coeff 1592 c$$$ if (rlatu(j).ge.0) then ! northern hemisphere 1593 c$$$ ! N.B: rlatu(:) contains decreasing values of latitude 1594 c$$$ ! starting from PI/2 to -PI/2 1595 c$$$ MONS_coeff=MONS_coeffN 1596 c$$$ else ! southern hemisphere 1597 c$$$ MONS_coeff=MONS_coeffS 1598 c$$$ endif 1599 c$$$ ! check if we should put subterranean ice 1600 c$$$ if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14% 1601 c$$$ ! compute depth at which ice lies: 1602 c$$$ val=MONS_d21(i,j)*MONS_coeff 1603 c$$$ ! compute val2= the diurnal skin depth of surface inertia 1604 c$$$ ! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1 1605 c$$$ val2=ith(i,j,1)*1.e-6*sqrt(88775./3.14159) 1606 c$$$ if (val.lt.val2) then 1607 c$$$ ! ice must be below the (surface inertia) diurnal skin depth 1608 c$$$ val=val2 1609 c$$$ endif 1610 c$$$ if (val.lt.layer(nsoilmx)) then ! subterranean ice 1611 c$$$ ! find the reference index iref that depth corresponds to 1612 c$$$ iref=0 1613 c$$$ do isoil=1,nsoilmx-1 1614 c$$$ if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1))) 1615 c$$$ & then 1616 c$$$ iref=isoil 1617 c$$$ exit 1618 c$$$ endif 1619 c$$$ enddo 1620 c$$$ ! Build "equivalent" thermal inertia for the mixed layer 1621 c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ 1622 c$$$ & (((val-layer(iref))/(ith(i,j,iref+1)**2))+ 1623 c$$$ & ((layer(iref+1)-val)/(iceith)**2))) 1624 c$$$ ! Set thermal inertia of lower layers 1625 c$$$ do isoil=iref+2,nsoilmx 1626 c$$$ ith(i,j,isoil)=iceith 1627 c$$$ enddo 1628 c$$$ endif ! of if (val.lt.layer(nsoilmx)) 1629 c$$$ endif ! of if (MONS_Hdn(i,j).lt.14.0) 1630 c$$$ enddo ! do j=1,jjp1 1631 c$$$ enddo ! do i=1,iip1 1632 c$$$ 1633 c$$$! Check: 1634 c$$$! do i=1,iip1 1635 c$$$! do j=1,jjp1 1636 c$$$! do isoil=1,nsoilmx 1637 c$$$! write(77,*) i,j,isoil," ",ith(i,j,isoil) 1638 c$$$! enddo 1639 c$$$! enddo 1640 c$$$! enddo 1641 c$$$ 1642 c$$$ ! recast ith() into ithfi() 1643 c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) 1644 c$$$ 1645 c$$$ else 1646 c$$$ write(*,*) ' Unknown (misspelled?) option!!!' 1647 end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ... 1266 else 1267 write(*,*) ' Unknown (misspelled?) option!!!' 1268 end if ! of if (trim(modif) .eq. '...') elseif ... 1648 1269 1649 1270 … … 1811 1432 end 1812 1433 1813 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1814 subroutine load_MONS_data(MONS_Hdn,MONS_d21)1815 use datafile_mod, only: datadir1816 implicit none1817 ! routine to load Benedicte Diez MONS dataset, fill in date in southern1818 ! polar region, and interpolate the result onto the GCM grid1819 #include"dimensions.h"1820 #include"paramet.h"1821 #include"comgeom.h"1822 ! arguments:1823 real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O1824 real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)1825 ! N.B MONS datasets should be of dimension (iip1,jjp1)1826 ! local variables:1827 character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"1828 character(len=88) :: txt ! to store some text1829 integer :: ierr,i,j1830 integer,parameter :: nblon=180 ! number of longitudes of MONS datasets1831 integer,parameter :: nblat=90 ! number of latitudes of MONS datasets1832 real :: pi1833 real :: longitudes(nblon) ! MONS dataset longitudes1834 real :: latitudes(nblat) ! MONS dataset latitudes1835 ! MONS dataset: mass fraction of H2O where H is assumed to be in H2O1836 real :: Hdn(nblon,nblat)1837 real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)1838 1839 ! Extended MONS dataset (for interp_horiz)1840 real :: Hdnx(nblon+1,nblat)1841 real :: d21x(nblon+1,nblat)1842 real :: lon_bound(nblon+1) ! longitude boundaries1843 real :: lat_bound(nblat-1) ! latitude boundaries1844 1845 ! 1. Initializations:1846 1847 write(*,*) "Loading MONS data"1848 1849 ! Open MONS datafile:1850 open(42,file=trim(datadir)//"/"//trim(filename),1851 & status="old",iostat=ierr)1852 if (ierr/=0) then1853 write(*,*) "Error in load_MONS_data:"1854 write(*,*) "Failed opening file ",1855 & trim(datadir)//"/"//trim(filename)1856 write(*,*)'Check that your path to datagcm:',trim(datadir)1857 write(*,*)' is correct. You can change it in callphys.def with:'1858 write(*,*)' datadir = /absolute/path/to/datagcm'1859 write(*,*)'If necessary ',trim(filename),1860 & ' (and other datafiles)'1861 write(*,*)' can be obtained online at:'1862 write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'1863 CALL ABORT1864 else ! skip first line of file (dummy read)1865 read(42,*) txt1866 endif1867 1868 pi=2.*asin(1.)1869 1870 !2. Load MONS data (on MONS grid)1871 do j=1,nblat1872 do i=1,nblon1873 ! swap latitude index so latitudes go from north pole to south pole:1874 read(42,*) latitudes(nblat-j+1),longitudes(i),1875 & Hdn(i,nblat-j+1),d21(i,nblat-j+1)1876 ! multiply d21 by 10 to convert from g/cm2 to kg/m21877 d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.01878 enddo1879 enddo1880 close(42)1881 1882 ! there is unfortunately no d21 data for latitudes -77 to -901883 ! so we build some by linear interpolation between values at -751884 ! and assuming d21=0 at the pole1885 do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-751886 do i=1,nblon1887 d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)1888 enddo1889 enddo1890 1891 ! 3. Build extended MONS dataset & boundaries (for interp_horiz)1892 ! longitude boundaries (in radians):1893 do i=1,nblon1894 ! NB: MONS data is every 2 degrees in longitude1895 lon_bound(i)=(longitudes(i)+1.0)*pi/180.01896 enddo1897 ! extra 'modulo' value1898 lon_bound(nblon+1)=lon_bound(1)+2.0*pi1899 1900 ! latitude boundaries (in radians):1901 do j=1,nblat-11902 ! NB: Mons data is every 2 degrees in latitude1903 lat_bound(j)=(latitudes(j)-1.0)*pi/180.01904 enddo1905 1906 ! MONS datasets:1907 do j=1,nblat1908 Hdnx(1:nblon,j)=Hdn(1:nblon,j)1909 Hdnx(nblon+1,j)=Hdnx(1,j)1910 d21x(1:nblon,j)=d21(1:nblon,j)1911 d21x(nblon+1,j)=d21x(1,j)1912 enddo1913 1914 ! Interpolate onto GCM grid1915 call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,1916 & lon_bound,lat_bound,rlonu,rlatv)1917 call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,1918 & lon_bound,lat_bound,rlonu,rlatv)1919 1920 end subroutine -
trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F
r305 r699 11 11 #include "dimensions.h" 12 12 #include "dimphys.h" 13 #include "comgeomfi.h"13 !#include "comgeomfi.h" 14 14 #include "surfdat.h" 15 15 #include "planete.h" … … 111 111 . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) 112 112 c 113 c Lecture des latitudes (coordonnees): 114 c 115 ierr = NF_INQ_VARID (nid, "latitude", nvarid) 116 IF (ierr.NE.NF_NOERR) THEN 117 PRINT*, 'phyetat0: Le champ <latitude> est absent' 118 CALL abort 119 ENDIF 120 #ifdef NC_DOUBLE 121 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lati) 122 #else 123 ierr = NF_GET_VAR_REAL(nid, nvarid, lati) 124 #endif 125 IF (ierr.NE.NF_NOERR) THEN 126 PRINT*, 'phyetat0: Lecture echouee pour <latitude>' 127 CALL abort 128 ENDIF 129 c 130 c Lecture des longitudes (coordonnees): 131 c 132 ierr = NF_INQ_VARID (nid, "longitude", nvarid) 133 IF (ierr.NE.NF_NOERR) THEN 134 PRINT*, 'phyetat0: Le champ <longitude> est absent' 135 CALL abort 136 ENDIF 137 #ifdef NC_DOUBLE 138 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, long) 139 #else 140 ierr = NF_GET_VAR_REAL(nid, nvarid, long) 141 #endif 142 IF (ierr.NE.NF_NOERR) THEN 143 PRINT*, 'phyetat0: Lecture echouee pour <longitude>' 144 CALL abort 145 ENDIF 146 c 147 c Lecture des aires des mailles: 148 c 149 ierr = NF_INQ_VARID (nid, "area", nvarid) 150 IF (ierr.NE.NF_NOERR) THEN 151 PRINT*, 'phyetat0: Le champ <area> est absent' 152 CALL abort 153 ENDIF 154 #ifdef NC_DOUBLE 155 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, area) 156 #else 157 ierr = NF_GET_VAR_REAL(nid, nvarid, area) 158 #endif 159 IF (ierr.NE.NF_NOERR) THEN 160 PRINT*, 'phyetat0: Lecture echouee pour <area>' 161 CALL abort 162 ENDIF 163 xmin = 1.0E+20 164 xmax = -1.0E+20 165 xmin = MINVAL(area) 166 xmax = MAXVAL(area) 167 PRINT*,'Aires des mailles <area>:', xmin, xmax 113 c Read latitudes (coordinates): No need, these are provided by the dynamics 114 c 115 ! ierr=nf90_inq_varid(nid,"latitude",nvarid) 116 ! IF (ierr.NE.nf90_noerr) THEN 117 ! PRINT*, 'phyetat0: Le champ <latitude> est absent' 118 ! write(*,*)trim(nf90_strerror(ierr)) 119 ! CALL abort 120 ! ENDIF 121 ! ierr=nf90_get_var(nid,nvarid,lati) 122 ! IF (ierr.NE.nf90_noerr) THEN 123 ! PRINT*, 'phyetat0: Lecture echouee pour <latitude>' 124 ! write(*,*)trim(nf90_strerror(ierr)) 125 ! CALL abort 126 ! ENDIF 127 c 128 c read longitudes (coordinates): No need, these are provided by the dynamics 129 c 130 ! ierr=nf90_inq_varid(nid,"longitude",nvarid) 131 ! IF (ierr.NE.nf90_noerr) THEN 132 ! PRINT*, 'phyetat0: Le champ <longitude> est absent' 133 ! write(*,*)trim(nf90_strerror(ierr)) 134 ! CALL abort 135 ! ENDIF 136 ! ierr=nf90_get_var(nid,nvarid,long) 137 ! IF (ierr.NE.nf90_noerr) THEN 138 ! PRINT*, 'phyetat0: Lecture echouee pour <longitude>' 139 ! write(*,*)trim(nf90_strerror(ierr)) 140 ! CALL abort 141 ! ENDIF 142 c 143 c Read areas of meshes: No need, these are provided by the dynamics 144 c 145 ! ierr=nf90_inq_varid(nid,"area",nvarid) 146 ! IF (ierr.NE.nf90_noerr) THEN 147 ! PRINT*, 'phyetat0: Le champ <area> est absent' 148 ! write(*,*)trim(nf90_strerror(ierr)) 149 ! CALL abort 150 ! ENDIF 151 ! ierr=nf90_get_var(nid,nvarid,area) 152 ! IF (ierr.NE.nf90_noerr) THEN 153 ! PRINT*, 'phyetat0: Lecture echouee pour <area>' 154 ! write(*,*)trim(nf90_strerror(ierr)) 155 ! CALL abort 156 ! ENDIF 157 ! xmin = 1.0E+20 158 ! xmax = -1.0E+20 159 ! xmin = MINVAL(area) 160 ! xmax = MAXVAL(area) 161 ! PRINT*,'Aires des mailles <area>:', xmin, xmax 168 162 c 169 163 c Lecture du geopotentiel au sol:
Note: See TracChangeset
for help on using the changeset viewer.
