Changeset 3198 for trunk/LMDZ.PLUTO


Ignore:
Timestamp:
Feb 2, 2024, 3:08:36 PM (10 months ago)
Author:
tbertrand
Message:

Cleaning newstart.F and adapting it to Pluto + small adjustments in the physics
TB

Location:
trunk/LMDZ.PLUTO/libf
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F

    r3195 r3198  
    22     &     date,tsurf,tsoil,emis,q2,
    33     &     t,ucov,vcov,ps,h,phisold_newgrid,
    4      &     q,qsurf,surfith,nid,
    5      &     rnat,pctsrf_sic,!tslab,tsea_ice,sea_ice,
    6      &     du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
     4     &     q,qsurf,surfith,nid)
     5!     &     rnat,pctsrf_sic) !tslab,tsea_ice,sea_ice,
     6! TB24     &     du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
    77
    88      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
    9       USE tracer_h, ONLY: igcm_haze
     9      USE tracer_h, ONLY: igcm_n2
    1010      USE infotrac, ONLY: tname, nqtot
    1111!      USE slab_ice_h, ONLY: noceanmx
    … …  
    7070      REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot)
    7171      ! REAL,INTENT(OUT) :: tslab(ngrid,nslay)
    72       REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid)
     72      !REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid)
    7373      ! REAL,INTENT(OUT) :: tsea_ice(ngrid),sea_ice(ngrid)
    7474c     REAL phisfi(ngrid)
    75       REAL,INTENT(OUT):: du_nonoro_gwd(ngrid,llm)
    76       REAL,INTENT(OUT):: dv_nonoro_gwd(ngrid,llm)
    77       REAL,INTENT(OUT):: east_gwstress(ngrid,llm)
    78       REAL,INTENT(OUT):: west_gwstress(ngrid,llm)
     75
     76c TB24
     77c      REAL,INTENT(OUT):: du_nonoro_gwd(ngrid,llm)
     78c      REAL,INTENT(OUT):: dv_nonoro_gwd(ngrid,llm)
     79c      REAL,INTENT(OUT):: east_gwstress(ngrid,llm)
     80c      REAL,INTENT(OUT):: west_gwstress(ngrid,llm)
    7981
    8082      INTEGER i,j,l
    … …  
    100102      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
    101103      ! real tslabS(iip1,jjp1,nslay)
    102       real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
    103       real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
    104       real du_nonoro_gwdS(iip1,jjp1,llm),dv_nonoro_gwdS(iip1,jjp1,llm)
    105       real east_gwstressS(iip1,jjp1,llm),west_gwstressS(iip1,jjp1,llm)
     104      !real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
     105      !real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
     106
     107!TB24
     108!      real du_nonoro_gwdS(iip1,jjp1,llm),dv_nonoro_gwdS(iip1,jjp1,llm)
     109!      real east_gwstressS(iip1,jjp1,llm),west_gwstressS(iip1,jjp1,llm)
    106110
    107111      real ptotal, n2icetotal
    … …  
    133137      real, dimension(:,:,:,:), allocatable :: qold
    134138      ! real, dimension(:,:,:), allocatable :: tslabold
    135       real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold
    136       real, dimension(:,:), allocatable :: tsea_iceold,sea_iceold
    137       real,allocatable :: du_nonoro_gwdold(:,:,:)
    138       real,allocatable :: dv_nonoro_gwdold(:,:,:)
    139       real,allocatable :: east_gwstressold(:,:,:)
    140       real,allocatable :: west_gwstressold(:,:,:)
     139      !real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold
     140      !real, dimension(:,:), allocatable :: tsea_iceold,sea_iceold
     141
     142      ! TB24
     143!      real,allocatable :: du_nonoro_gwdold(:,:,:)
     144!      real,allocatable :: dv_nonoro_gwdold(:,:,:)
     145!      real,allocatable :: east_gwstressold(:,:,:)
     146!      real,allocatable :: west_gwstressold(:,:,:)
    141147
    142148      real tab_cntrl(100)
    … …  
    289295      allocate(qsurfold(imold+1,jmold+1,nqtot))
    290296      ! allocate(tslabold(imold+1,jmold+1,nslay))
    291       allocate(rnatold(imold+1,jmold+1))
    292       allocate(pctsrf_sicold(imold+1,jmold+1))
    293       allocate(tsea_iceold(imold+1,jmold+1))
    294       allocate(sea_iceold(imold+1,jmold+1))
    295 
    296       allocate(du_nonoro_gwdold(imold+1,jmold+1,lmold))
    297       allocate(dv_nonoro_gwdold(imold+1,jmold+1,lmold))
    298       allocate(east_gwstressold(imold+1,jmold+1,lmold))
    299       allocate(west_gwstressold(imold+1,jmold+1,lmold))
     297      !allocate(rnatold(imold+1,jmold+1))
     298      !allocate(pctsrf_sicold(imold+1,jmold+1))
     299      !allocate(tsea_iceold(imold+1,jmold+1))
     300      !allocate(sea_iceold(imold+1,jmold+1))
     301
     302      !TB24
     303!      allocate(du_nonoro_gwdold(imold+1,jmold+1,lmold))
     304!      allocate(dv_nonoro_gwdold(imold+1,jmold+1,lmold))
     305!      allocate(east_gwstressold(imold+1,jmold+1,lmold))
     306!      allocate(west_gwstressold(imold+1,jmold+1,lmold))
    300307
    301308      allocate(var (imold+1,jmold+1,llm))
    … …  
    960967      ENDDO ! of DO iq=1,nqtot
    961968
    962 ! Non-orographic GWs:
    963       write(*,*)"lect_start_archive: loading du_nonoro_gwd"
    964       ierr = NF_INQ_VARID (nid,"du_nonoro_gwd", nvarid)
    965       IF (ierr .NE. NF_NOERR) THEN
    966          PRINT*, "lect_start_archive: Field <du_nonoro_gwd> not found"
    967          PRINT*, "Setting it to zero"
    968          du_nonoro_gwdold(:,:,:)=0
    969       ELSE
    970 #ifdef NC_DOUBLE
    971         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,du_nonoro_gwdold)
    972 #else
    973         ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,du_nonoro_gwdold)
    974 #endif
    975         IF (ierr .NE. NF_NOERR) THEN
    976           PRINT*, "lect_start_archive: Failed loading <du_nonoro_gwd>"
    977           CALL abort
    978         ENDIF
    979       ENDIF
    980 
    981       write(*,*)"lect_start_archive: loading dv_nonoro_gwd"
    982       ierr = NF_INQ_VARID (nid,"dv_nonoro_gwd", nvarid)
    983       IF (ierr .NE. NF_NOERR) THEN
    984          PRINT*, "lect_start_archive: Field <dv_nonoro_gwd> not found"
    985          PRINT*, "Setting it to zero"
    986          dv_nonoro_gwdold(:,:,:)=0
    987       ELSE
    988 #ifdef NC_DOUBLE
    989         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,dv_nonoro_gwdold)
    990 #else
    991         ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,dv_nonoro_gwdold)
    992 #endif
    993         IF (ierr .NE. NF_NOERR) THEN
    994           PRINT*, "lect_start_archive: Failed loading <dv_nonoro_gwd>"
    995           CALL abort
    996         ENDIF
    997       ENDIF
    998 
    999       write(*,*)"lect_start_archive: loading east_gwstress"
    1000       ierr = NF_INQ_VARID (nid,"east_gwstress", nvarid)
    1001       IF (ierr .NE. NF_NOERR) THEN
    1002          PRINT*, "lect_start_archive: Field <east_gwstress> not found"
    1003          PRINT*, "Setting it to zero"
    1004          east_gwstressold(:,:,:)=0
    1005       ELSE
    1006 #ifdef NC_DOUBLE
    1007         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,east_gwstressold)
    1008 #else
    1009         ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,east_gwstressold)
    1010 #endif
    1011         IF (ierr .NE. NF_NOERR) THEN
    1012           PRINT*, "lect_start_archive: Failed loading <east_gwstress>"
    1013           CALL abort
    1014         ENDIF
    1015       ENDIF
    1016 
    1017       write(*,*)"lect_start_archive: loading west_gwstress"
    1018       ierr = NF_INQ_VARID (nid,"west_gwstress", nvarid)
    1019       IF (ierr .NE. NF_NOERR) THEN
    1020          PRINT*, "lect_start_archive: Field <west_gwstress> not found"
    1021          PRINT*, "Setting it to zero"
    1022          west_gwstressold(:,:,:)=0
    1023       ELSE
    1024 #ifdef NC_DOUBLE
    1025         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,west_gwstressold)
    1026 #else
    1027         ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,west_gwstressold)
    1028 #endif
    1029         IF (ierr .NE. NF_NOERR) THEN
    1030           PRINT*, "lect_start_archive: Failed loading <west_gwstress>"
    1031           CALL abort
    1032         ENDIF
    1033       ENDIF
     969! Non-orographic GWs: TB24 rm
     970
     971!      write(*,*)"lect_start_archive: loading du_nonoro_gwd"
     972!      ierr = NF_INQ_VARID (nid,"du_nonoro_gwd", nvarid)
     973!      IF (ierr .NE. NF_NOERR) THEN
     974!         PRINT*, "lect_start_archive: Field <du_nonoro_gwd> not found"
     975!         PRINT*, "Setting it to zero"
     976!         du_nonoro_gwdold(:,:,:)=0
     977!      ELSE
     978!#ifdef NC_DOUBLE
     979!        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,du_nonoro_gwdold)
     980!#else
     981!        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,du_nonoro_gwdold)
     982!#endif
     983!        IF (ierr .NE. NF_NOERR) THEN
     984!          PRINT*, "lect_start_archive: Failed loading <du_nonoro_gwd>"
     985!          CALL abort
     986!        ENDIF
     987!      ENDIF
     988
     989!      write(*,*)"lect_start_archive: loading dv_nonoro_gwd"
     990!      ierr = NF_INQ_VARID (nid,"dv_nonoro_gwd", nvarid)
     991!      IF (ierr .NE. NF_NOERR) THEN
     992!         PRINT*, "lect_start_archive: Field <dv_nonoro_gwd> not found"
     993!         PRINT*, "Setting it to zero"
     994!         dv_nonoro_gwdold(:,:,:)=0
     995!      ELSE
     996!#ifdef NC_DOUBLE
     997!        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,dv_nonoro_gwdold)
     998!#else
     999!        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,dv_nonoro_gwdold)
     1000!#endif
     1001!        IF (ierr .NE. NF_NOERR) THEN
     1002!          PRINT*, "lect_start_archive: Failed loading <dv_nonoro_gwd>"
     1003!          CALL abort
     1004!        ENDIF
     1005!      ENDIF
     1006
     1007!      write(*,*)"lect_start_archive: loading east_gwstress"
     1008!      ierr = NF_INQ_VARID (nid,"east_gwstress", nvarid)
     1009!      IF (ierr .NE. NF_NOERR) THEN
     1010!         PRINT*, "lect_start_archive: Field <east_gwstress> not found"
     1011!         PRINT*, "Setting it to zero"
     1012!         east_gwstressold(:,:,:)=0
     1013!      ELSE
     1014!#ifdef NC_DOUBLE
     1015!        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,east_gwstressold)
     1016!#else
     1017!        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,east_gwstressold)
     1018!#endif
     1019!        IF (ierr .NE. NF_NOERR) THEN
     1020!          PRINT*, "lect_start_archive: Failed loading <east_gwstress>"
     1021!          CALL abort
     1022!        ENDIF
     1023!      ENDIF
     1024
     1025!      write(*,*)"lect_start_archive: loading west_gwstress"
     1026!      ierr = NF_INQ_VARID (nid,"west_gwstress", nvarid)
     1027!      IF (ierr .NE. NF_NOERR) THEN
     1028!         PRINT*, "lect_start_archive: Field <west_gwstress> not found"
     1029!         PRINT*, "Setting it to zero"
     1030!         west_gwstressold(:,:,:)=0
     1031!      ELSE
     1032!#ifdef NC_DOUBLE
     1033!        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,west_gwstressold)
     1034!#else
     1035!!        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,west_gwstressold)
     1036!#endif
     1037!        IF (ierr .NE. NF_NOERR) THEN
     1038!          PRINT*, "lect_start_archive: Failed loading <west_gwstress>"
     1039!          CALL abort
     1040!        ENDIF
     1041!      ENDIF
    10341042
    10351043!=======================================================================
    … …  
    10521060     &                   rlonuold,rlatvold,rlonu,rlatv)
    10531061
    1054 ! N2 ice is now in qsurf(igcm_haze)
     1062! N2 ice is now in qsurf(igcm_n2)
    10551063!      call interp_horiz (n2iceold,n2ices,imold,jmold,iim,jjm,1,
    10561064!     &                   rlonuold,rlatvold,rlonu,rlatv)
    … …  
    10971105      END DO
    10981106      n2icetotal = 0.
    1099       if (igcm_haze.ne.0) then
     1107      if (igcm_n2.ne.0) then
    11001108        ! recast surface N2 ice on new grid
    1101         call interp_horiz(qsurfold(1,1,igcm_haze),
    1102      &                  qsurfs(1,1,igcm_haze),
     1109        call interp_horiz(qsurfold(1,1,igcm_n2),
     1110     &                  qsurfs(1,1,igcm_n2),
    11031111     &                  imold,jmold,iim,jjm,1,
    11041112     &                  rlonuold,rlatvold,rlonu,rlatv)
    … …  
    11061114         DO i=1,iim
    11071115           !n2icetotal = n2icetotal + n2iceS(i,j)*aire(i,j)
    1108            n2icetotal=n2icetotal+qsurfS(i,j,igcm_haze)*aire(i,j)
     1116           n2icetotal=n2icetotal+qsurfS(i,j,igcm_n2)*aire(i,j)
    11091117         END DO
    11101118       END DO
    … …  
    13151323      write (*,*) 'lect_start_archive: t ', t(1,jjp1,1)  ! INFO
    13161324
    1317 ! Non-orographic GW
    1318       call interp_vert
    1319      &    (du_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
    1320      &     psold,(imold+1)*(jmold+1))
    1321       call interp_horiz(var,du_nonoro_gwdS,imold,jmold,iim,jjm,llm,
    1322      &                   rlonuold,rlatvold,rlonu,rlatv)
    1323       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd)
    1324 
    1325       call interp_vert
    1326      &    (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
    1327      &     psold,(imold+1)*(jmold+1))
    1328       call interp_horiz(var,dv_nonoro_gwdS,imold,jmold,iim,jjm,llm,
    1329      &                   rlonuold,rlatvold,rlonu,rlatv)
    1330       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd)
    1331 
    1332       call interp_vert
    1333      &    (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
    1334      &     psold,(imold+1)*(jmold+1))
    1335       call interp_horiz(var,east_gwstressS,imold,jmold,iim,jjm,llm,
    1336      &                   rlonuold,rlatvold,rlonu,rlatv)
    1337       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress)
    1338 
    1339       call interp_vert
    1340      &    (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
    1341      &     psold,(imold+1)*(jmold+1))
    1342       call interp_horiz(var,west_gwstressS,imold,jmold,iim,jjm,llm,
    1343      &                   rlonuold,rlatvold,rlonu,rlatv)
    1344       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,west_gwstressS,west_gwstress)
    1345 
     1325! Non-orographic GW TB24 rm
     1326!      call interp_vert
     1327!     &    (du_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
     1328!     &     psold,(imold+1)*(jmold+1))
     1329!      call interp_horiz(var,du_nonoro_gwdS,imold,jmold,iim,jjm,llm,
     1330!     &                   rlonuold,rlatvold,rlonu,rlatv)
     1331!      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd)
     1332
     1333!      call interp_vert
     1334!     &    (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
     1335!     &     psold,(imold+1)*(jmold+1))
     1336!      call interp_horiz(var,dv_nonoro_gwdS,imold,jmold,iim,jjm,llm,
     1337!     &                   rlonuold,rlatvold,rlonu,rlatv)
     1338!      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd)
     1339!
     1340!      call interp_vert
     1341!     &    (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
     1342!     &     psold,(imold+1)*(jmold+1))
     1343!      call interp_horiz(var,east_gwstressS,imold,jmold,iim,jjm,llm,
     1344!     &                   rlonuold,rlatvold,rlonu,rlatv)
     1345!      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress)
     1346!
     1347!      call interp_vert
     1348!     &    (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
     1349!     &     psold,(imold+1)*(jmold+1))
     1350!      call interp_horiz(var,west_gwstressS,imold,jmold,iim,jjm,llm,
     1351!     &                   rlonuold,rlatvold,rlonu,rlatv)
     1352!      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,west_gwstressS,west_gwstress)
     1353!
    13461354c q2 : pbl wind variance
    13471355      write (*,*) 'lect_start_archive: q2old ', q2old (1,2,1)  ! INFO
    … …  
    14541462
    14551463!      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,n2ices,n2ice)
    1456 ! no need to transfer "n2ice" any more; it is in qsurf(igcm_haze)
     1464! no need to transfer "n2ice" any more; it is in qsurf(igcm_n2)
    14571465
    14581466      endif !! if nqtot .ne. 0
    … …  
    15011509      deallocate(var,varp1)
    15021510      ! deallocate(tslabold)
    1503       deallocate(rnatold)
    1504       deallocate(pctsrf_sicold)
     1511      !deallocate(rnatold)
     1512      !deallocate(pctsrf_sicold)
    15051513      ! deallocate(tsea_iceold)
    15061514      ! deallocate(sea_iceold)

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

    r3184 r3198  
    1919     &                              is_master
    2020      use infotrac, only: infotrac_init, nqtot, tname
    21       USE tracer_h, ONLY: igcm_n2_ice, igcm_h2o_vap, igcm_h2o_ice
     21      USE tracer_h, ONLY: igcm_n2
    2222      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
    2323      USE surfdat_h, ONLY: phisfi, albedodat,
    2424     &                     zmea, zstd, zsig, zgam, zthe
    25       USE nonoro_gwd_ran_mod, ONLY: du_nonoro_gwd, dv_nonoro_gwd,
    26      &                              east_gwstress, west_gwstress
     25! TB24      USE nonoro_gwd_ran_mod, ONLY: du_nonoro_gwd, dv_nonoro_gwd,
     26!     &                              east_gwstress, west_gwstress
    2727      use datafile_mod, only: datadir, surfdir
    2828      use ioipsl_getin_p_mod, only: getin_p
    … …  
    3131      use iostart, only: open_startphy
    3232!      use slab_ice_h, only:noceanmx
    33       USE ocean_slab_mod, ONLY: nslay
     33!      USE ocean_slab_mod, ONLY: nslay
    3434      use filtreg_mod, only: inifilr
    3535      USE mod_const_mpi, ONLY: COMM_LMDZ
    … …  
    180180!      real :: time_step,t_ops,t_wrt
    181181!      CHARACTER*80 :: visu_file
     182      !TB24
     183      CALL conf_gcm( 99, .TRUE. )
     184
    182185
    183186      cpp    = 0.
    … …  
    208211      allocate(tsoil(ngridmx,nsoilmx))
    209212      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
    210       allocate(tslab(ngridmx,nsoilmx))
    211213     
    212214c=======================================================================
    … …  
    351353        CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
    352354     .        nqtot,day_ini,time,
    353      .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
    354      .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
    355      .        sea_ice)
     355     .        tsurf,tsoil,emis,q2,qsurf)    !) ! temporary modif by RDW
     356!     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
     357!     .        sea_ice)
    356358
    357359        ! copy albedo and soil thermal inertia on (local) physics grid
    … …  
    545547     &   date,tsurf,tsoil,emis,q2,
    546548     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
    547      &   surfith,nid,
    548      &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice,
    549      &   du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
     549     &   surfith,nid)
     550
     551!     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     552! TB24     &   du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
    550553        write(*,*) "OK, read start_archive file"
    551554        ! copy soil thermal inertia
    … …  
    576579      write(*,*) 'flat : no topography ("aquaplanet")'
    577580      write(*,*) 'set_ps_to_preff : used if changing preff with topo'
    578       write(*,*) 'nuketharsis : no Tharsis bulge'
    579       write(*,*) 'bilball : uniform albedo and thermal inertia'
    580       write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
    581581      write(*,*) 'qname : change tracer name'
    582582      write(*,*) 't=profile  : read temperature profile in profile.in'
    … …  
    585585      write(*,*) 'qs=x : give a uniform value to a surface tracer'
    586586      write(*,*) 'q=profile    : specify a profile for a tracer'
    587 !      write(*,*) 'ini_q : tracers initialisation for chemistry, water an
    588 !     $d ice   '
    589 !      write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and
    590 !     $ice '
    591 !      write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on
    592 !     $ly '
    593       write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45'
    594       write(*,*) 'watercapn : H20 ice on permanent N polar cap '
    595       write(*,*) 'watercaps : H20 ice on permanent S polar cap '
    596       write(*,*) 'noacglac  : H2O ice across Noachis Terra'
    597       write(*,*) 'oborealis : H2O ice across Vastitas Borealis'
    598       write(*,*) 'iceball   : Thick ice layer all over surface'
    599       write(*,*) 'supercontinent: Create a continent of given Ab and TI'
    600       write(*,*) 'wetstart  : start with a wet atmosphere'
    601       write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
    602       write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
    603      $ profile (lat-alt) and winds set to zero'
    604       write(*,*) 'coldstart : Start X K above the N2 frost point and
    605      $set wind to zero (assumes 100% N2)'
    606       write(*,*) 'n2ice=0 : remove N2 polar cap'
    607       write(*,*) 'ptot : change total pressure'
    608       write(*,*) 'emis : change surface emissivity'
    609       write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
    610      &face values'
    611       write(*,*) 'slab_ocean_0 : initialisation of slab
    612      $ocean variables'
    613       write(*,*) 'chemistry_ini : initialisation of chemical profiles'
    614587
    615588        write(*,*)
    … …  
    673646           enddo
    674647          enddo
    675 
    676 c       'nuketharsis : no tharsis bulge for Early Mars'
    677 c       ---------------------------------------------
    678         else if (trim(modif) .eq. 'nuketharsis') then
    679 
    680            DO j=1,jjp1       
    681               DO i=1,iim
    682                  ig=1+(j-2)*iim +i
    683                  if(j.eq.1) ig=1
    684                  if(j.eq.jjp1) ig=ngridmx
    685 
    686                  fact1=(((rlonv(i)*180./pi)+100)**2 +
    687      &                (rlatu(j)*180./pi)**2)/65**2
    688                  fact2=exp( -fact1**2.5 )
    689 
    690                  phis(i,j) = phis(i,j) - (phis(i,j)+4000.*g)*fact2
    691 
    692 !                 if(phis(i,j).gt.2500.*g)then
    693 !                    if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
    694 !                       phis(i,j)=2500.*g
    695 !                    endif
    696 !                 endif
    697 
    698               ENDDO
    699            ENDDO
    700           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    701 
    702 
    703 c       bilball : uniform albedo, thermal inertia
    704 c       -----------------------------------------
    705         else if (trim(modif) .eq. 'bilball') then
    706           write(*,*) 'constante albedo and iner.therm:'
    707           write(*,*) 'New value for albedo (ex: 0.25) ?'
    708  101      read(*,*,iostat=ierr) alb_bb
    709           if(ierr.ne.0) goto 101
    710           write(*,*)
    711           write(*,*) ' uniform albedo (new value):',alb_bb
    712           write(*,*)
    713 
    714           write(*,*) 'New value for thermal inertia (eg: 247) ?'
    715  102      read(*,*,iostat=ierr) ith_bb
    716           if(ierr.ne.0) goto 102
    717           write(*,*) 'uniform thermal inertia (new value):',ith_bb
    718           DO j=1,jjp1
    719              DO i=1,iip1
    720                 alb(i,j) = alb_bb        ! albedo
    721                 do isoil=1,nsoilmx
    722                   ith(i,j,isoil) = ith_bb        ! thermal inertia
    723                 enddo
    724              END DO
    725           END DO
    726 !          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
    727           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    728           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    729 
    730 c       coldspole : sous-sol de la calotte sud toujours froid
    731 c       -----------------------------------------------------
    732         else if (trim(modif) .eq. 'coldspole') then
    733           write(*,*)'new value for the subsurface temperature',
    734      &   ' beneath the permanent southern polar cap ? (eg: 141 K)'
    735  103      read(*,*,iostat=ierr) tsud
    736           if(ierr.ne.0) goto 103
    737           write(*,*)
    738           write(*,*) ' new value of the subsurface temperature:',tsud
    739 c         nouvelle temperature sous la calotte permanente
    740           do l=2,nsoilmx
    741                tsoil(ngridmx,l) =  tsud
    742           end do
    743 
    744 
    745           write(*,*)'new value for the albedo',
    746      &   'of the permanent southern polar cap ? (eg: 0.75)'
    747  104      read(*,*,iostat=ierr) albsud
    748           if(ierr.ne.0) goto 104
    749           write(*,*)
    750 
    751 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    752 c         Option 1:  only the albedo of the pole is modified :   
    753           albfi(ngridmx)=albsud
    754           write(*,*) 'ig=',ngridmx,'   albedo perennial cap ',
    755      &    albfi(ngridmx)
    756 
    757 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    758 c          Option 2 A haute resolution : coordonnee de la vrai calotte ~   
    759 c           DO j=1,jjp1
    760 c             DO i=1,iip1
    761 c                ig=1+(j-2)*iim +i
    762 c                if(j.eq.1) ig=1
    763 c                if(j.eq.jjp1) ig=ngridmx
    764 c                if ((rlatu(j)*180./pi.lt.-84.).and.
    765 c     &            (rlatu(j)*180./pi.gt.-91.).and.
    766 c     &            (rlonv(i)*180./pi.gt.-91.).and.
    767 c     &            (rlonv(i)*180./pi.lt.0.))         then
    768 cc    albedo de la calotte permanente fixe a albsud
    769 c                   alb(i,j)=albsud
    770 c                   write(*,*) 'lat=',rlatu(j)*180./pi,
    771 c     &                      ' lon=',rlonv(i)*180./pi
    772 cc     fin de la condition sur les limites de la calotte permanente
    773 c                end if
    774 c             ENDDO
    775 c          ENDDO
    776 c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    777 
    778 c         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    779 
    780 
    781 c       ptot : Modification of the total pressure: ice + current atmosphere
    782 c       -------------------------------------------------------------------
    783         else if (trim(modif).eq.'ptot') then
    784 
    785           ! check if we have a n2_ice surface tracer:
    786           if (igcm_n2_ice.eq.0) then
    787             write(*,*) " No surface N2 ice !"
    788             write(*,*) " only atmospheric pressure will be considered!"
    789           endif
    790 c         calcul de la pression totale glace + atm actuelle
    791           patm=0.
    792           airetot=0.
    793           pcap=0.
    794           DO j=1,jjp1
    795              DO i=1,iim
    796                 ig=1+(j-2)*iim +i
    797                 if(j.eq.1) ig=1
    798                 if(j.eq.jjp1) ig=ngridmx
    799                 patm = patm + ps(i,j)*aire(i,j)
    800                 airetot= airetot + aire(i,j)
    801                 if (igcm_n2_ice.ne.0) then
    802                   !pcap = pcap + aire(i,j)*n2ice(ig)*g
    803                   pcap = pcap + aire(i,j)*qsurf(ig,igcm_n2_ice)*g
    804                 endif
    805              ENDDO
    806           ENDDO
    807           ptoto = pcap + patm
    808 
    809           print*,'Current total pressure at surface (n2 ice + atm) ',
    810      &       ptoto/airetot
    811 
    812           print*,'new value?'
    813           read(*,*) ptotn
    814           ptotn=ptotn*airetot
    815           patmn=ptotn-pcap
    816           print*,'ptoto,patm,ptotn,patmn'
    817           print*,ptoto,patm,ptotn,patmn
    818           print*,'Mult. factor for pressure (atm only)', patmn/patm
    819           do j=1,jjp1
    820              do i=1,iip1
    821                 ps(i,j)=ps(i,j)*patmn/patm
    822              enddo
    823           enddo
    824 
    825 
    826 
    827 c        Correction pour la conservation des traceurs
    828          yes=' '
    829          do while ((yes.ne.'y').and.(yes.ne.'n'))
    830             write(*,*) 'Do you wish to conserve tracer total mass (y)',
    831      &              ' or tracer mixing ratio (n) ?'
    832              read(*,fmt='(a)') yes
    833          end do
    834 
    835          if (yes.eq.'y') then
    836            write(*,*) 'OK : conservation of tracer total mass'
    837            DO iq =1, nqtot
    838              DO l=1,llm
    839                DO j=1,jjp1
    840                   DO i=1,iip1
    841                     q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
    842                   ENDDO
    843                ENDDO
    844              ENDDO
    845            ENDDO
    846           else
    847             write(*,*) 'OK : conservation of tracer mixing ratio'
    848           end if
    849 
    850 c        Correction pour la pression au niveau de la mer
    851          yes=' '
    852          do while ((yes.ne.'y').and.(yes.ne.'n'))
    853             write(*,*) 'Do you wish fix pressure at sea level (y)',
    854      &              ' or not (n) ?'
    855              read(*,fmt='(a)') yes
    856          end do
    857 
    858          if (yes.eq.'y') then
    859            write(*,*) 'Value?'
    860                 read(*,*,iostat=ierr) psea
    861              DO i=1,iip1
    862                DO j=1,jjp1
    863                     ps(i,j)=psea
    864 
    865                   ENDDO
    866                ENDDO
    867                 write(*,*) 'psea=',psea
    868           else
    869             write(*,*) 'no'
    870           end if
    871 
    872 
    873 c       emis : change surface emissivity (added by RW)
    874 c       ----------------------------------------------
    875         else if (trim(modif).eq.'emis') then
    876 
    877            print*,'new value?'
    878            read(*,*) emisread
    879 
    880            do i=1,ngridmx
    881               emis(i)=emisread
    882            enddo
    883648
    884649c       qname : change tracer name
    … …  
    1057822               write(*,*)'No modifications to tracer ',trim(tname(iq))
    1058823             endif
     824
     825
     826        endif
    1059827             
    1060 
    1061 c      wetstart : wet atmosphere with a north to south gradient
    1062 c      --------------------------------------------------------
    1063        else if (trim(modif) .eq. 'wetstart') then
    1064         ! check that there is indeed a water vapor tracer
    1065         if (igcm_h2o_vap.eq.0) then
    1066           write(*,*) "No water vapour tracer! Can't use this option"
    1067           stop
    1068         endif
    1069           DO l=1,llm
    1070             DO j=1,jjp1
    1071               DO i=1,iip1
    1072                 q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
    1073               ENDDO
    1074             ENDDO
    1075           ENDDO
    1076 
    1077          write(*,*) 'Water mass mixing ratio at north pole='
    1078      *               ,q(1,1,1,igcm_h2o_vap)
    1079          write(*,*) '---------------------------south pole='
    1080      *               ,q(1,jjp1,1,igcm_h2o_vap)
    1081 
    1082 c      noglacier : remove tropical water ice (to initialize high res sim)
    1083 c      --------------------------------------------------
    1084         else if (trim(modif) .eq. 'noglacier') then
    1085            if (igcm_h2o_ice.eq.0) then
    1086              write(*,*) "No water ice tracer! Can't use this option"
    1087              stop
    1088            endif
    1089            do ig=1,ngridmx
    1090              j=(ig-2)/iim +2
    1091               if(ig.eq.1) j=1
    1092               write(*,*) 'OK: remove surface ice for |lat|<45'
    1093               if (abs(rlatu(j)*180./pi).lt.45.) then
    1094                    qsurf(ig,igcm_h2o_ice)=0.
    1095               end if
    1096            end do
    1097 
    1098 
    1099 c      watercapn : H20 ice on permanent northern cap
    1100 c      --------------------------------------------------
    1101         else if (trim(modif) .eq. 'watercapn') then
    1102            if (igcm_h2o_ice.eq.0) then
    1103              write(*,*) "No water ice tracer! Can't use this option"
    1104              stop
    1105            endif
    1106 
    1107            DO j=1,jjp1       
    1108               DO i=1,iim
    1109                  ig=1+(j-2)*iim +i
    1110                  if(j.eq.1) ig=1
    1111                  if(j.eq.jjp1) ig=ngridmx
    1112 
    1113                  if (rlatu(j)*180./pi.gt.80.) then
    1114                     qsurf(ig,igcm_h2o_ice)=3.4e3
    1115                     !do isoil=1,nsoilmx
    1116                     !   ith(i,j,isoil) = 18000. ! thermal inertia
    1117                     !enddo
    1118                    write(*,*)'     ==> Ice mesh North boundary (deg)= ',
    1119      &                   rlatv(j-1)*180./pi
    1120                  end if
    1121               ENDDO
    1122            ENDDO
    1123            CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1124 
    1125 c      watercaps : H20 ice on permanent southern cap
    1126 c      -------------------------------------------------
    1127         else if (trim(modif) .eq. 'watercaps') then
    1128            if (igcm_h2o_ice.eq.0) then
    1129               write(*,*) "No water ice tracer! Can't use this option"
    1130               stop
    1131            endif
    1132 
    1133            DO j=1,jjp1       
    1134               DO i=1,iim
    1135                  ig=1+(j-2)*iim +i
    1136                  if(j.eq.1) ig=1
    1137                  if(j.eq.jjp1) ig=ngridmx
    1138 
    1139                  if (rlatu(j)*180./pi.lt.-80.) then
    1140                     qsurf(ig,igcm_h2o_ice)=3.4e3
    1141                     !do isoil=1,nsoilmx
    1142                     !   ith(i,j,isoil) = 18000. ! thermal inertia
    1143                     !enddo
    1144                    write(*,*)'     ==> Ice mesh South boundary (deg)= ',
    1145      &                   rlatv(j-1)*180./pi
    1146                  end if
    1147               ENDDO
    1148            ENDDO
    1149            CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1150 
    1151 
    1152 c       noacglac : H2O ice across highest terrain
    1153 c       --------------------------------------------
    1154         else if (trim(modif) .eq. 'noacglac') then
    1155            if (igcm_h2o_ice.eq.0) then
    1156              write(*,*) "No water ice tracer! Can't use this option"
    1157              stop
    1158            endif
    1159           DO j=1,jjp1       
    1160              DO i=1,iim
    1161                 ig=1+(j-2)*iim +i
    1162                 if(j.eq.1) ig=1
    1163                 if(j.eq.jjp1) ig=ngridmx
    1164 
    1165                 if(phis(i,j).gt.1000.*g)then
    1166                     alb(i,j) = 0.5 ! snow value
    1167                     do isoil=1,nsoilmx
    1168                        ith(i,j,isoil) = 18000. ! thermal inertia
    1169                        ! this leads to rnat set to 'ocean' in physiq.F90
    1170                        ! actually no, because it is soil not surface
    1171                     enddo
    1172                 endif
    1173              ENDDO
    1174           ENDDO
    1175           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1176           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1177           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1178 
    1179 
    1180 
    1181 c       oborealis : H2O oceans across Vastitas Borealis
    1182 c       -----------------------------------------------
    1183         else if (trim(modif) .eq. 'oborealis') then
    1184            if (igcm_h2o_ice.eq.0) then
    1185              write(*,*) "No water ice tracer! Can't use this option"
    1186              stop
    1187            endif
    1188           DO j=1,jjp1       
    1189              DO i=1,iim
    1190                 ig=1+(j-2)*iim +i
    1191                 if(j.eq.1) ig=1
    1192                 if(j.eq.jjp1) ig=ngridmx
    1193 
    1194                 if(phis(i,j).lt.-4000.*g)then
    1195 !                if( (phis(i,j).lt.-4000.*g)
    1196 !     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
    1197 !                    phis(i,j)=-8000.0*g ! proper ocean
    1198                     phis(i,j)=-4000.0*g ! scanty ocean
    1199 
    1200                     alb(i,j) = 0.07 ! oceanic value
    1201                     do isoil=1,nsoilmx
    1202                        ith(i,j,isoil) = 18000. ! thermal inertia
    1203                        ! this leads to rnat set to 'ocean' in physiq.F90
    1204                        ! actually no, because it is soil not surface
    1205                     enddo
    1206                 endif
    1207              ENDDO
    1208           ENDDO
    1209           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1210           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1211           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1212 
    1213 c       iborealis : H2O ice in Northern plains
    1214 c       --------------------------------------
    1215         else if (trim(modif) .eq. 'iborealis') then
    1216            if (igcm_h2o_ice.eq.0) then
    1217              write(*,*) "No water ice tracer! Can't use this option"
    1218              stop
    1219            endif
    1220           DO j=1,jjp1       
    1221              DO i=1,iim
    1222                 ig=1+(j-2)*iim +i
    1223                 if(j.eq.1) ig=1
    1224                 if(j.eq.jjp1) ig=ngridmx
    1225 
    1226                 if(phis(i,j).lt.-4000.*g)then
    1227                    !qsurf(ig,igcm_h2o_ice)=1.e3
    1228                    qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
    1229                 endif
    1230              ENDDO
    1231           ENDDO
    1232           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1233           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1234           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1235 
    1236 
    1237 c       oceanball : H2O liquid everywhere
    1238 c       ----------------------------
    1239         else if (trim(modif) .eq. 'oceanball') then
    1240            if (igcm_h2o_ice.eq.0) then
    1241              write(*,*) "No water ice tracer! Can't use this option"
    1242              stop
    1243            endif
    1244           DO j=1,jjp1       
    1245              DO i=1,iim
    1246                 ig=1+(j-2)*iim +i
    1247                 if(j.eq.1) ig=1
    1248                 if(j.eq.jjp1) ig=ngridmx
    1249 
    1250                 qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
    1251                 qsurf(ig,igcm_h2o_ice)=0.0
    1252                 alb(i,j) = 0.07 ! ocean value
    1253 
    1254                 do isoil=1,nsoilmx
    1255                    ith(i,j,isoil) = 18000. ! thermal inertia
    1256                    !ith(i,j,isoil) = 50. ! extremely low for test
    1257                    ! this leads to rnat set to 'ocean' in physiq.F90
    1258                 enddo
    1259 
    1260              ENDDO
    1261           ENDDO
    1262           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1263           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1264           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    1265 
    1266 c       iceball : H2O ice everywhere
    1267 c       ----------------------------
    1268         else if (trim(modif) .eq. 'iceball') then
    1269            if (igcm_h2o_ice.eq.0) then
    1270              write(*,*) "No water ice tracer! Can't use this option"
    1271              stop
    1272            endif
    1273           DO j=1,jjp1       
    1274              DO i=1,iim
    1275                 ig=1+(j-2)*iim +i
    1276                 if(j.eq.1) ig=1
    1277                 if(j.eq.jjp1) ig=ngridmx
    1278 
    1279                 qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
    1280                 qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
    1281                 alb(i,j) = 0.6 ! ice albedo value
    1282 
    1283                 do isoil=1,nsoilmx
    1284                    !ith(i,j,isoil) = 18000. ! thermal inertia
    1285                    ! this leads to rnat set to 'ocean' in physiq.F90
    1286                 enddo
    1287 
    1288              ENDDO
    1289           ENDDO
    1290           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1291           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1292 
    1293 c       supercontinent : H2O ice everywhere
    1294 c       ----------------------------
    1295         else if (trim(modif) .eq. 'supercontinent') then
    1296           write(*,*) 'Minimum longitude (-180,180)'
    1297           read(*,*) val
    1298           write(*,*) 'Maximum longitude (-180,180)'
    1299           read(*,*) val2
    1300           write(*,*) 'Minimum latitude (-90,90)'
    1301           read(*,*) val3
    1302           write(*,*) 'Maximum latitude (-90,90)'
    1303           read(*,*) val4
    1304 
    1305           do j=1,jjp1
    1306             do i=1,iip1
    1307               ig=1+(j-2)*iim +i
    1308               if(j.eq.1) ig=1
    1309               if(j.eq.jjp1) ig=ngridmx
    1310 
    1311 c             Supercontinent:
    1312               if (((rlatu(j)*180./pi.le.val4).and.
    1313      &            (rlatu(j)*180./pi.ge.val3).and.
    1314      &            (rlonv(i)*180./pi.le.val2).and.
    1315      &            (rlonv(i)*180./pi.ge.val))) then
    1316 
    1317                 rnat(ig)=1.
    1318                 alb(i,j) = 0.3
    1319                 do isoil=1,nsoilmx
    1320                   ith(i,j,isoil) = 2000.
    1321                 enddo
    1322 c             Ocean:
    1323               else
    1324                 rnat(ig)=0.
    1325                 alb(i,j) = 0.07
    1326                 do isoil=1,nsoilmx
    1327                   ith(i,j,isoil) = 18000.
    1328                 enddo
    1329               end if
    1330 
    1331             enddo
    1332           enddo
    1333           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1334           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    1335 
    1336 c       isotherm : Isothermal temperatures and no winds
    1337 c       -----------------------------------------------
    1338         else if (trim(modif) .eq. 'isotherm') then
    1339 
    1340           write(*,*)'Isothermal temperature of the atmosphere,
    1341      &           surface and subsurface'
    1342           write(*,*) 'Value of this temperature ? :'
    1343  203      read(*,*,iostat=ierr) Tiso
    1344           if(ierr.ne.0) goto 203
    1345 
    1346           tsurf(1:ngridmx)=Tiso
    1347          
    1348           tsoil(1:ngridmx,1:nsoilmx)=Tiso
    1349          
    1350           Tset(1:iip1,1:jjp1,1:llm)=Tiso
    1351           flagtset=.true.
    1352 
    1353           t(1:iip1,1:jjp1,1:llm)=Tiso
    1354           !! otherwise hydrost. integrations below
    1355           !! use the wrong temperature
    1356           !! -- NB: Tset might be useless
    1357        
    1358           ucov(1:iip1,1:jjp1,1:llm)=0
    1359           vcov(1:iip1,1:jjm,1:llm)=0
    1360           q2(1:ngridmx,1:llm+1)=0
    1361 
    1362 c       radequi : Radiative equilibrium profile of temperatures and no winds
    1363 c       --------------------------------------------------------------------
    1364         else if (trim(modif) .eq. 'radequi') then
    1365 
    1366           write(*,*)'radiative equilibrium temperature profile'       
    1367 
    1368           do ig=1, ngridmx
    1369              teque= 335.0-60.0*sin(latitude(ig))*sin(latitude(ig))-
    1370      &            10.0*cos(latitude(ig))*cos(latitude(ig))
    1371              tsurf(ig) = MAX(220.0,teque)
    1372           end do
    1373           do l=2,nsoilmx
    1374              do ig=1, ngridmx
    1375                 tsoil(ig,l) = tsurf(ig)
    1376              end do
    1377           end do
    1378           DO j=1,jjp1
    1379              DO i=1,iim
    1380                 Do l=1,llm
    1381                    teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))-
    1382      &                  10.0*cos(rlatu(j))*cos(rlatu(j))
    1383                    Tset(i,j,l)=MAX(220.0,teque)
    1384                 end do
    1385              end do
    1386           end do
    1387           flagtset=.true.
    1388           ucov(1:iip1,1:jjp1,1:llm)=0
    1389           vcov(1:iip1,1:jjm,1:llm)=0
    1390           q2(1:ngridmx,1:llm+1)=0
    1391 
    1392 c       coldstart : T set 1K above N2 frost point and no winds
    1393 c       ------------------------------------------------
    1394         else if (trim(modif) .eq. 'coldstart') then
    1395 
    1396           write(*,*)'set temperature of the atmosphere,'
    1397      &,'surface and subsurface how many degrees above N2 frost point?'
    1398  204      read(*,*,iostat=ierr) Tabove
    1399           if(ierr.ne.0) goto 204
    1400 
    1401             DO j=1,jjp1
    1402              DO i=1,iim
    1403                 ig=1+(j-2)*iim +i
    1404                 if(j.eq.1) ig=1
    1405                 if(j.eq.jjp1) ig=ngridmx
    1406             tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
    1407              END DO
    1408             END DO
    1409           do l=1,nsoilmx
    1410             do ig=1, ngridmx
    1411               tsoil(ig,l) = tsurf(ig)
    1412             end do
    1413           end do
    1414           DO j=1,jjp1
    1415            DO i=1,iim
    1416             Do l=1,llm
    1417                pp = aps(l) +bps(l)*ps(i,j)
    1418                Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
    1419             end do
    1420            end do
    1421           end do
    1422 
    1423           flagtset=.true.
    1424           ucov(1:iip1,1:jjp1,1:llm)=0
    1425           vcov(1:iip1,1:jjm,1:llm)=0
    1426           q2(1:ngridmx,1:llm+1)=0
    1427 
    1428 
    1429 c       n2ice=0 : remove N2 polar ice caps'
    1430 c       ------------------------------------------------
    1431         else if (trim(modif) .eq. 'n2ice=0') then
    1432          ! check that there is indeed a n2_ice tracer ...
    1433           if (igcm_n2_ice.ne.0) then
    1434            do ig=1,ngridmx
    1435               !n2ice(ig)=0
    1436               qsurf(ig,igcm_n2_ice)=0
    1437               emis(ig)=emis(ngridmx/2)
    1438            end do
    1439           else
    1440             write(*,*) "Can't remove N2 ice!! (no n2_ice tracer)"
    1441           endif
    1442        
    1443 !       therm_ini_s: (re)-set soil thermal inertia to reference surface values
    1444 !       ----------------------------------------------------------------------
    1445 
    1446         else if (trim(modif) .eq. 'therm_ini_s') then
    1447 !          write(*,*)"surfithfi(1):",surfithfi(1)
    1448           do isoil=1,nsoilmx
    1449             inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
    1450           enddo
    1451           write(*,*)'OK: Soil thermal inertia has been reset to referenc
    1452      &e surface values'
    1453 !          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
    1454           ithfi(:,:)=inertiedat(:,:)
    1455          ! recast ithfi() onto ith()
    1456          call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    1457 ! Check:
    1458 !         do i=1,iip1
    1459 !           do j=1,jjp1
    1460 !             do isoil=1,nsoilmx
    1461 !               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
    1462 !             enddo
    1463 !           enddo
    1464 !         enddo
    1465 
    1466 
    1467 
    1468 c       slab_ocean_initialisation
    1469 c       ------------------------------------------------
    1470         else if (trim(modif) .eq. 'slab_ocean_0') then
    1471         write(*,*)'OK: initialisation of slab ocean'
    1472 
    1473       DO ig=1, ngridmx
    1474          rnat(ig)=1.
    1475          tslab(ig,1)=0.
    1476          tslab(ig,2)=0.
    1477          tsea_ice(ig)=0.
    1478          sea_ice(ig)=0.
    1479          pctsrf_sic(ig)=0.
    1480          
    1481          if(ithfi(ig,1).GT.10000.)then
    1482            rnat(ig)=0.
    1483            phisfi(ig)=0.
    1484            tsea_ice(ig)=273.15-1.8
    1485            tslab(ig,1)=tsurf(ig)
    1486            tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
    1487           endif
    1488 
    1489       ENDDO
    1490           CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
    1491 
    1492 
    1493 c       chemistry_initialisation
    1494 c       ------------------------------------------------
    1495         else if (trim(modif) .eq. 'chemistry_ini') then
    1496         write(*,*)'OK: initialisation of chemical profiles'
    1497 
    1498 
    1499           call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
    1500      &                          1, 0)
    1501 
    1502          ! We want to have the very same value at lon -180 and lon 180
    1503           do l = 1,llm
    1504              do j = 1,jjp1
    1505                 do iq = 1,nqtot
    1506                    q(iip1,j,l,iq) = q(1,j,l,iq)
    1507                 end do
    1508              end do
    1509           end do
    1510 
    1511 
    1512         else
    1513           write(*,*) '       Unknown (misspelled?) option!!!'
    1514         end if ! of if (trim(modif) .eq. '...') elseif ...
    1515 
    1516 
    1517 
    1518828       enddo ! of do ! infinite loop on liste of changes
    1519829
    … …  
    1524834c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
    1525835c=======================================================================
    1526       DO ig=1, ngridmx
    1527          totalfrac(ig)=0.5
    1528          DO l=1,llm
    1529             cloudfrac(ig,l)=0.5
    1530          ENDDO
     836!      DO ig=1, ngridmx
     837!         totalfrac(ig)=0.5
     838!         DO l=1,llm
     839!            cloudfrac(ig,l)=0.5
     840!         ENDDO
    1531841! Ehouarn, march 2012: also add some initialisation for hice
    1532          hice(ig)=0.0
    1533       ENDDO
     842!         hice(ig)=0.0
     843!      ENDDO
    1534844     
    1535845c========================================================
    … …  
    1651961      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
    1652962     &                dtphys,real(day_ini),
    1653      &                tsurf,tsoil,emis,q2,qsurf,
    1654      &                cloudfrac,totalfrac,hice,
    1655      &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     963     &                tsurf,tsoil,emis,q2,qsurf)
     964!     &                cloudfrac,totalfrac,hice,
     965!     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
    1656966
    1657967c=======================================================================

  • trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r3197 r3198  
    920920               print*,'and the surface albedo is taken equal to the first visible spectral value'
    921921
     922               albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
    922923               fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
     924               ! TB24:
     925               fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
    923926               fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
    924927               fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
    … …  
    16621665      if(callrad)then
    16631666
    1664          !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
     1667         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    16651668         !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
    16661669         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)

  • trunk/LMDZ.PLUTO/libf/phypluto/surfini.F

    r3195 r3198  
    4444      ! Step 2 : We get the bare ground albedo from the start files.
    4545      DO ig=1,ngrid
    46          albedo_bareground(ig)=albedodat(ig)
     46         albedo_bareground(ig)=0.1 ! TB24 albedodat(ig)
    4747         DO nw=1,L_NSPECTV
    48             albedo(ig,nw)=albedo_bareground(ig)
     48            albedo(ig,nw)=0.1 !albedo_bareground(ig)
    4949         ENDDO
    5050      ENDDO

  • trunk/LMDZ.PLUTO/libf/phypluto/turbdiff_mod.F90

    r3197 r3198  
    138138
    139139      IF (firstcall) THEN
     140         ivap=1 !TB24
     141         iliq=0
     142         iliq_surf=0
     143         iice_surf=0 ! simply to make the code legible                   
    140144         sensibFlux(:)=0.
    141145         firstcall=.false.
Note: See TracChangeset for help on using the changeset viewer.