Ignore:
Timestamp:
Jan 11, 2017, 3:33:51 PM (8 years ago)
Author:
jvatant
Message:

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F

    r1644 r1647  
    1919     &                              is_master
    2020      use infotrac, only: infotrac_init, nqtot, tname
    21       USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
    2221      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
    2322      USE surfdat_h, ONLY: phisfi, albedodat,
    … …  
    2827      use phyredem, only: physdem0, physdem1
    2928      use iostart, only: open_startphy
    30       use slab_ice_h, only:noceanmx
    3129      use filtreg_mod, only: inifilr
    3230      USE mod_const_mpi, ONLY: COMM_LMDZ
    … …  
    9593      REAL tsurf(ngridmx)        ! surface temperature
    9694      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
    97 !      REAL co2ice(ngridmx)        ! CO2 ice layer
    9895      REAL emis(ngridmx)        ! surface emissivity
    9996      real emisread             ! added by RW
    … …  
    109106      real mugaz ! molar mass of the atmosphere
    110107
    111       integer ierr
    112 
    113       REAL rnat(ngridmx)
    114       REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
    115       REAL pctsrf_sic(ngridmx) ! sea ice cover
    116       REAL tsea_ice(ngridmx) ! temperature sea_ice
    117       REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2)
     108      integer ierr
    118109
    119110c Variables on the new grid along scalar points
    … …  
    153144      logical :: flagtset=.false. ,  flagps0=.false.
    154145      real val, val2, val3, val4 ! to store temporary variables
    155       real :: iceith=2000 ! thermal inertia of subterranean ice
    156146
    157147      INTEGER :: itau
    … …  
    165155!     added by BC for equilibrium temperature startup
    166156      real teque
    167 
    168 !     added by BC for cloud fraction setup
    169       REAL hice(ngridmx),cloudfrac(ngridmx,llm)
    170       REAL totalfrac(ngridmx)
    171 
    172157
    173158!     added by RW for nuketharsis
    … …  
    209194      allocate(tsoil(ngridmx,nsoilmx))
    210195      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
    211       allocate(tslab(ngridmx,nsoilmx))
    212196     
    213197c=======================================================================
    … …  
    355339        CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
    356340     .        nqtot,day_ini,time,
    357      .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
    358      .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
    359      .        sea_ice)
     341     .        tsurf,tsoil,emis,q2,qsurf)
    360342
    361343        ! copy albedo and soil thermal inertia on (local) physics grid
    … …  
    549531     &   date,tsurf,tsoil,emis,q2,
    550532     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
    551      &   surfith,nid,
    552      &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     533     &   surfith,nid)
    553534        write(*,*) "OK, read start_archive file"
    554535        ! copy soil thermal inertia
    … …  
    784765        else if (trim(modif).eq.'ptot') then
    785766
    786           ! check if we have a co2_ice surface tracer:
    787           if (igcm_co2_ice.eq.0) then
    788             write(*,*) " No surface CO2 ice !"
    789             write(*,*) " only atmospheric pressure will be considered!"
    790           endif
    791767c         calcul de la pression totale glace + atm actuelle
    792768          patm=0.
    … …  
    800776                patm = patm + ps(i,j)*aire(i,j)
    801777                airetot= airetot + aire(i,j)
    802                 if (igcm_co2_ice.ne.0) then
    803                   !pcap = pcap + aire(i,j)*co2ice(ig)*g
    804                   pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2_ice)*g
    805                 endif
    806778             ENDDO
    807779          ENDDO
    808780          ptoto = pcap + patm
    809781
    810           print*,'Current total pressure at surface (co2 ice + atm) ',
     782          print*,'Current total pressure at surface (atm) ',
    811783     &       ptoto/airetot
    812784
    … …  
    10481020       else if (trim(modif) .eq. 'wetstart') then
    10491021        ! check that there is indeed a water vapor tracer
    1050         if (igcm_h2o_vap.eq.0) then
     1022
    10511023          write(*,*) "No water vapour tracer! Can't use this option"
    10521024          stop
    1053         endif
    1054           DO l=1,llm
    1055             DO j=1,jjp1
    1056               DO i=1,iip1
    1057                 q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
    1058               ENDDO
    1059             ENDDO
    1060           ENDDO
    1061 
    1062          write(*,*) 'Water mass mixing ratio at north pole='
    1063      *               ,q(1,1,1,igcm_h2o_vap)
    1064          write(*,*) '---------------------------south pole='
    1065      *               ,q(1,jjp1,1,igcm_h2o_vap)
    10661025
    10671026c      noglacier : remove tropical water ice (to initialize high res sim)
    10681027c      --------------------------------------------------
    10691028        else if (trim(modif) .eq. 'noglacier') then
    1070            if (igcm_h2o_ice.eq.0) then
     1029       
    10711030             write(*,*) "No water ice tracer! Can't use this option"
    10721031             stop
    1073            endif
    1074            do ig=1,ngridmx
    1075              j=(ig-2)/iim +2
    1076               if(ig.eq.1) j=1
    1077               write(*,*) 'OK: remove surface ice for |lat|<45'
    1078               if (abs(rlatu(j)*180./pi).lt.45.) then
    1079                    qsurf(ig,igcm_h2o_ice)=0.
    1080               end if
    1081            end do
     1032
    10821033
    10831034
    … …  
    10851036c      --------------------------------------------------
    10861037        else if (trim(modif) .eq. 'watercapn') then
    1087            if (igcm_h2o_ice.eq.0) then
     1038       
    10881039             write(*,*) "No water ice tracer! Can't use this option"
    10891040             stop
    1090            endif
    1091 
    1092            DO j=1,jjp1       
    1093               DO i=1,iim
    1094                  ig=1+(j-2)*iim +i
    1095                  if(j.eq.1) ig=1
    1096                  if(j.eq.jjp1) ig=ngridmx
    1097 
    1098                  if (rlatu(j)*180./pi.gt.80.) then
    1099                     qsurf(ig,igcm_h2o_ice)=3.4e3
    1100                     !do isoil=1,nsoilmx
    1101                     !   ith(i,j,isoil) = 18000. ! thermal inertia
    1102                     !enddo
    1103                    write(*,*)'     ==> Ice mesh North boundary (deg)= ',
    1104      &                   rlatv(j-1)*180./pi
    1105                  end if
    1106               ENDDO
    1107            ENDDO
    1108            CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1109 
    1110 c$$$           do ig=1,ngridmx
    1111 c$$$             j=(ig-2)/iim +2
    1112 c$$$              if(ig.eq.1) j=1
    1113 c$$$              if (rlatu(j)*180./pi.gt.80.) then
    1114 c$$$
    1115 c$$$                   qsurf(ig,igcm_h2o_ice)=1.e5
    1116 c$$$                   qsurf(ig,igcm_h2o_vap)=0.0!1.e5
    1117 c$$$
    1118 c$$$                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
    1119 c$$$     &              qsurf(ig,igcm_h2o_ice)
    1120 c$$$
    1121 c$$$                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
    1122 c$$$     &              rlatv(j)*180./pi
    1123 c$$$                end if
    1124 c$$$           enddo
    11251041
    11261042c      watercaps : H20 ice on permanent southern cap
    11271043c      -------------------------------------------------
    11281044        else if (trim(modif) .eq. 'watercaps') then
    1129            if (igcm_h2o_ice.eq.0) then
     1045
    11301046              write(*,*) "No water ice tracer! Can't use this option"
    11311047              stop
    1132            endif
    1133 
    1134            DO j=1,jjp1       
    1135               DO i=1,iim
    1136                  ig=1+(j-2)*iim +i
    1137                  if(j.eq.1) ig=1
    1138                  if(j.eq.jjp1) ig=ngridmx
    1139 
    1140                  if (rlatu(j)*180./pi.lt.-80.) then
    1141                     qsurf(ig,igcm_h2o_ice)=3.4e3
    1142                     !do isoil=1,nsoilmx
    1143                     !   ith(i,j,isoil) = 18000. ! thermal inertia
    1144                     !enddo
    1145                    write(*,*)'     ==> Ice mesh South boundary (deg)= ',
    1146      &                   rlatv(j-1)*180./pi
    1147                  end if
    1148               ENDDO
    1149            ENDDO
    1150            CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1151 
    1152 c$$$           do ig=1,ngridmx
    1153 c$$$               j=(ig-2)/iim +2
    1154 c$$$               if(ig.eq.1) j=1
    1155 c$$$               if (rlatu(j)*180./pi.lt.-80.) then
    1156 c$$$                  qsurf(ig,igcm_h2o_ice)=1.e5
    1157 c$$$                  qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
    1158 c$$$
    1159 c$$$                  write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
    1160 c$$$     &                 qsurf(ig,igcm_h2o_ice)
    1161 c$$$                  write(*,*)'     ==> Ice mesh North boundary (deg)= ',
    1162 c$$$     &                 rlatv(j-1)*180./pi
    1163 c$$$               end if
    1164 c$$$           enddo
    1165 
    11661048
    11671049c       noacglac : H2O ice across highest terrain
    11681050c       --------------------------------------------
    11691051        else if (trim(modif) .eq. 'noacglac') then
    1170            if (igcm_h2o_ice.eq.0) then
     1052
    11711053             write(*,*) "No water ice tracer! Can't use this option"
    11721054             stop
    1173            endif
    1174           DO j=1,jjp1       
    1175              DO i=1,iim
    1176                 ig=1+(j-2)*iim +i
    1177                 if(j.eq.1) ig=1
    1178                 if(j.eq.jjp1) ig=ngridmx
    1179 
    1180                 if(phis(i,j).gt.1000.*g)then
    1181                     alb(i,j) = 0.5 ! snow value
    1182                     do isoil=1,nsoilmx
    1183                        ith(i,j,isoil) = 18000. ! thermal inertia
    1184                        ! this leads to rnat set to 'ocean' in physiq.F90
    1185                        ! actually no, because it is soil not surface
    1186                     enddo
    1187                 endif
    1188              ENDDO
    1189           ENDDO
    1190           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1191           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1192           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1193 
    1194 
    11951055
    11961056c       oborealis : H2O oceans across Vastitas Borealis
    11971057c       -----------------------------------------------
    11981058        else if (trim(modif) .eq. 'oborealis') then
    1199            if (igcm_h2o_ice.eq.0) then
     1059
    12001060             write(*,*) "No water ice tracer! Can't use this option"
    12011061             stop
    1202            endif
    1203           DO j=1,jjp1       
    1204              DO i=1,iim
    1205                 ig=1+(j-2)*iim +i
    1206                 if(j.eq.1) ig=1
    1207                 if(j.eq.jjp1) ig=ngridmx
    1208 
    1209                 if(phis(i,j).lt.-4000.*g)then
    1210 !                if( (phis(i,j).lt.-4000.*g)
    1211 !     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
    1212 !                    phis(i,j)=-8000.0*g ! proper ocean
    1213                     phis(i,j)=-4000.0*g ! scanty ocean
    1214 
    1215                     alb(i,j) = 0.07 ! oceanic value
    1216                     do isoil=1,nsoilmx
    1217                        ith(i,j,isoil) = 18000. ! thermal inertia
    1218                        ! this leads to rnat set to 'ocean' in physiq.F90
    1219                        ! actually no, because it is soil not surface
    1220                     enddo
    1221                 endif
    1222              ENDDO
    1223           ENDDO
    1224           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1225           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1226           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1227 
     1062             
    12281063c       iborealis : H2O ice in Northern plains
    12291064c       --------------------------------------
    12301065        else if (trim(modif) .eq. 'iborealis') then
    1231            if (igcm_h2o_ice.eq.0) then
     1066
    12321067             write(*,*) "No water ice tracer! Can't use this option"
    12331068             stop
    1234            endif
    1235           DO j=1,jjp1       
    1236              DO i=1,iim
    1237                 ig=1+(j-2)*iim +i
    1238                 if(j.eq.1) ig=1
    1239                 if(j.eq.jjp1) ig=ngridmx
    1240 
    1241                 if(phis(i,j).lt.-4000.*g)then
    1242                    !qsurf(ig,igcm_h2o_ice)=1.e3
    1243                    qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
    1244                 endif
    1245              ENDDO
    1246           ENDDO
    1247           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1248           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1249           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1250 
    12511069
    12521070c       oceanball : H2O liquid everywhere
    12531071c       ----------------------------
    12541072        else if (trim(modif) .eq. 'oceanball') then
    1255            if (igcm_h2o_ice.eq.0) then
     1073
    12561074             write(*,*) "No water ice tracer! Can't use this option"
    12571075             stop
    1258            endif
    1259           DO j=1,jjp1       
    1260              DO i=1,iim
    1261                 ig=1+(j-2)*iim +i
    1262                 if(j.eq.1) ig=1
    1263                 if(j.eq.jjp1) ig=ngridmx
    1264 
    1265                 qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
    1266                 qsurf(ig,igcm_h2o_ice)=0.0
    1267                 alb(i,j) = 0.07 ! ocean value
    1268 
    1269                 do isoil=1,nsoilmx
    1270                    ith(i,j,isoil) = 18000. ! thermal inertia
    1271                    !ith(i,j,isoil) = 50. ! extremely low for test
    1272                    ! this leads to rnat set to 'ocean' in physiq.F90
    1273                 enddo
    1274 
    1275              ENDDO
    1276           ENDDO
    1277           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1278           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1279           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    12801076
    12811077c       iceball : H2O ice everywhere
    12821078c       ----------------------------
    12831079        else if (trim(modif) .eq. 'iceball') then
    1284            if (igcm_h2o_ice.eq.0) then
     1080
    12851081             write(*,*) "No water ice tracer! Can't use this option"
    12861082             stop
    1287            endif
    1288           DO j=1,jjp1       
    1289              DO i=1,iim
    1290                 ig=1+(j-2)*iim +i
    1291                 if(j.eq.1) ig=1
    1292                 if(j.eq.jjp1) ig=ngridmx
    1293 
    1294                 qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
    1295                 qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
    1296                 alb(i,j) = 0.6 ! ice albedo value
    1297 
    1298                 do isoil=1,nsoilmx
    1299                    !ith(i,j,isoil) = 18000. ! thermal inertia
    1300                    ! this leads to rnat set to 'ocean' in physiq.F90
    1301                 enddo
    1302 
    1303              ENDDO
    1304           ENDDO
    1305           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1306           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    13071083
    13081084c       supercontinent : H2O ice everywhere
    13091085c       ----------------------------
    13101086        else if (trim(modif) .eq. 'supercontinent') then
    1311           write(*,*) 'Minimum longitude (-180,180)'
    1312           read(*,*) val
    1313           write(*,*) 'Maximum longitude (-180,180)'
    1314           read(*,*) val2
    1315           write(*,*) 'Minimum latitude (-90,90)'
    1316           read(*,*) val3
    1317           write(*,*) 'Maximum latitude (-90,90)'
    1318           read(*,*) val4
    1319 
    1320           do j=1,jjp1
    1321             do i=1,iip1
    1322               ig=1+(j-2)*iim +i
    1323               if(j.eq.1) ig=1
    1324               if(j.eq.jjp1) ig=ngridmx
    1325 
    1326 c             Supercontinent:
    1327               if (((rlatu(j)*180./pi.le.val4).and.
    1328      &            (rlatu(j)*180./pi.ge.val3).and.
    1329      &            (rlonv(i)*180./pi.le.val2).and.
    1330      &            (rlonv(i)*180./pi.ge.val))) then
    1331 
    1332                 rnat(ig)=1.
    1333                 alb(i,j) = 0.3
    1334                 do isoil=1,nsoilmx
    1335                   ith(i,j,isoil) = 2000.
    1336                 enddo
    1337 c             Ocean:
    1338               else
    1339                 rnat(ig)=0.
    1340                 alb(i,j) = 0.07
    1341                 do isoil=1,nsoilmx
    1342                   ith(i,j,isoil) = 18000.
    1343                 enddo
    1344               end if
    1345 
    1346             enddo
    1347           enddo
    1348           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1349           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
     1087 
     1088             write(*,*) "No water ice tracer! Can't use this option"
     1089             stop
    13501090
    13511091c       isotherm : Isothermal temperatures and no winds
    … …  
    14451185c       ------------------------------------------------
    14461186        else if (trim(modif) .eq. 'co2ice=0') then
    1447          ! check that there is indeed a co2_ice tracer ...
    1448           if (igcm_co2_ice.ne.0) then
    1449            do ig=1,ngridmx
    1450               !co2ice(ig)=0
    1451               qsurf(ig,igcm_co2_ice)=0
    1452               emis(ig)=emis(ngridmx/2)
    1453            end do
    1454           else
    14551187            write(*,*) "Can't remove CO2 ice!! (no co2_ice tracer)"
    1456           endif
     1188       
    14571189       
    14581190!       therm_ini_s: (re)-set soil thermal inertia to reference surface values
    … …  
    14841216c       ------------------------------------------------
    14851217        else if (trim(modif) .eq. 'slab_ocean_0') then
    1486         write(*,*)'OK: initialisation of slab ocean'
    1487 
    1488       DO ig=1, ngridmx
    1489          rnat(ig)=1.
    1490          tslab(ig,1)=0.
    1491          tslab(ig,2)=0.
    1492          tsea_ice(ig)=0.
    1493          sea_ice(ig)=0.
    1494          pctsrf_sic(ig)=0.
    1495          
    1496          if(ithfi(ig,1).GT.10000.)then
    1497            rnat(ig)=0.
    1498            phisfi(ig)=0.
    1499            tsea_ice(ig)=273.15-1.8
    1500            tslab(ig,1)=tsurf(ig)
    1501            tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
    1502           endif
    1503 
    1504       ENDDO
    1505           CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
    1506 
    1507 
     1218       
     1219           write(*,*) "No ocean yet on Titan ! Can't use this option"
     1220           stop
    15081221
    15091222        else
    … …  
    15171230 999  continue
    15181231
    1519  
    1520 c=======================================================================
    1521 c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
    1522 c=======================================================================
    1523       DO ig=1, ngridmx
    1524          totalfrac(ig)=0.5
    1525          DO l=1,llm
    1526             cloudfrac(ig,l)=0.5
    1527          ENDDO
    1528 ! Ehouarn, march 2012: also add some initialisation for hice
    1529          hice(ig)=0.0
    1530       ENDDO
    15311232     
    1532 c========================================================
    1533 
    1534 !      DO ig=1,ngridmx
    1535 !         IF(tsurf(ig) .LT. 273.16-1.8) THEN
    1536 !            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
    1537 !            hice(ig)=min(hice(ig),1.0)
    1538 !         ENDIF
    1539 !      ENDDO
    1540 
    1541 
    1542 
    15431233
    15441234c=======================================================================
    … …  
    16471337      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
    16481338     &                dtphys,real(day_ini),
    1649      &                tsurf,tsoil,emis,q2,qsurf,
    1650      &                cloudfrac,totalfrac,hice,
    1651      &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     1339     &                tsurf,tsoil,emis,q2,qsurf)
    16521340
    16531341c=======================================================================
Note: See TracChangeset for help on using the changeset viewer.