Changeset 1974 for trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars

Timestamp:

Jul 18, 2018, 4:48:34 PM (6 years ago)

Author:

mvals

Message:

Mars GCM:
Integration of the detached dust layer parametrizations (rocket dust storm, slope wind lifting, CW, and dust injection scheme, DB).
Still experimental, default behaviour (rdstorm=.false., dustinjection=0) identical to previous revision.
NB: Updated newstart requires an updated "surface.nc" containing the "hmons" field.
EM+MV

Location:

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars

Files:

• avg_horiz_mod.F (added)
• datareadnc.F (modified) (6 diffs)
• mvc_horiz_mod.F (added)
• newstart.F (modified) (5 diffs)

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F

@@ -1 +1 @@
 c=======================================================================
       SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
-     &                    zmea,zstd,zsig,zgam,zthe)
+     &                    zmea,zstd,zsig,zgam,zthe,
+     &                    hmons,summit,zavg)
 c=======================================================================
 c
@@ -42 +43 @@
 c=======================================================================
 
-      use ioipsl_getincom, only: getin
-      USE comconst_mod, ONLY: g,pi
+! to use  'getin'
+      use ioipsl_getincom, only: getin 
+      use comconst_mod, only: g,pi
       use datafile_mod, only: datadir
+      use avg_horiz_mod, only: avg_horiz
+      use mvc_horiz_mod, only: mvc_horiz
 
       implicit none
@@ -74 +78 @@
       REAL,intent(out) :: zgam(imdp1*jmdp1)
       REAL,intent(out) :: zthe(imdp1*jmdp1)
-
+      REAL,intent(out) :: hmons(imdp1*jmdp1) !CW17,361*180 hmons
+      REAL,intent(out) :: summit(imdp1*jmdp1) !CW17,361*180 hmons
+      REAL,intent(out) :: zavg(imdp1*jmdp1) !CW17,361*180 hmons
+      
 ! Local variables:
       REAL        zdata(imd*jmdp1)
@@ -95 +102 @@
 
       CHARACTER*20 string
-      DIMENSION string(0:4)
+      DIMENSION string(0:6)
 
 
@@ -233 +240 @@
           call grid_noro1(360, 180, longitude, latitude, zdata,
      .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
+     
+          !CW17
+          call avg_horiz(imd,jmdp1,iim,jjm,longitude,
+     .                latitude,rlonv,rlatu,zdata,pfield)
+
+          do j=1,jjp1 !CW17 49, iimp1=65, the last column = first column
+             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
+          enddo 
+          do i=1,iimp1*jjp1
+c             if (pfield(i) .ne. -999999.) then
+                zavg(i) = pfield(i)
+c             else
+c               zavg(i)=-999999.
+c             endif
+          enddo
 
       endif
@@ -297 +319 @@
 c    FIN
 c-----------------------------------------------------------------------
-
+ 
+c=======================================================================
+c<<<<<<<<<<<<<<<<<<<<<<<add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! name of dataset
+      string(5) = 'hmons' !subgrid hmons
+c     the following data could be useful in future, but currently,
+c     we don't need to put them into restartfi.nc
+      string(6) = 'summit' !subgrid summit
+c     string(7) = 'base'   !base of summit (not subgrid)
+c     string(8) = 'bottom' !subgrid base
+      do k=5,6
+          write(*,*) 'string',k,string(k)  
+c-----------------------------------------------------------------------
+c    Lecture NetCDF  
+c-----------------------------------------------------------------------
+
+          ierr = NF_INQ_VARID (unit, string(k), nvarid)
+          if (ierr.ne.nf_noerr) then
+            write(*,*) 'datareadnc error, cannot find ',trim(string(k))
+            write(*,*) ' in file ',trim(datadir),'/surface.nc'
+            stop
+          endif
+
+c-----------------------------------------------------------------------
+c    initialisation
+c-----------------------------------------------------------------------
+          call initial0(iimp1*jjp1,pfield)
+          call initial0(imd*jmdp1,zdata)
+          call initial0(imdp1*jmdp1,zdataS)
+
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
+#else
+          ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
+#endif
+          if (ierr.ne.nf_noerr) then
+            write(*,*) 'datareadnc: error failed loading ',
+     .                  trim(string(k))
+            stop
+          endif
+
+c-----------------------------------------------------------------------
+c   Passage de zdata en grille (imdp1 jmdp1)
+c-----------------------------------------------------------------------
+          do j=1,jmdp1     !  180
+              do i=1,imd   !  360
+                  !copy zdata to zdataS, line by line
+                  !      i+ 361 *(j-1)          i+  360*(j-1)
+                  zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
+              enddo
+              ! the last column = the first column of zdata
+              zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
+          enddo
+
+          if (k .eq. 5 .or. k .eq. 6) then 
+              ! hmons, summit, keep the maximum of subgrids 
+              call mvc_horiz(imd,jmdp1,iim,jjm,longitude,
+     .                latitude,rlonv,rlatu,zdata,pfield)
+          endif
+
+
+c-----------------------------------------------------------------------
+c    Periodicite    
+c-----------------------------------------------------------------------
+
+          do j=1,jjp1 ! 49, iimp1=65, the last column = first column
+             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
+          enddo 
+ 
+c-----------------------------------------------------------------------
+c    Sauvegarde des champs    
+c-----------------------------------------------------------------------
+
+          if (k.eq.5) then ! hmons
+             do i=1,iimp1*jjp1
+                  if (pfield(i) .ne. -999999.) then
+                    hmons(i) = pfield(i)
+                  else
+                    hmons(i)=0.
+                  endif
+             enddo
+          endif
+          
+          if (k.eq.6) then  ! summit
+             do i=1,iimp1*jjp1
+                    summit(i) = pfield(i)
+             enddo
+          endif
+
+      enddo
+c<<<<<<<<<<<<<<<<<<<<<<<<<done add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>
+c=======================================================================
       END

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F

@@ -27 +27 @@
       use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
      &                     albedodat, z0_default, qsurf, tsurf,
-     &                     co2ice, emis
+     &                     co2ice, emis, hmons
       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
       use control_mod, only: day_step, iphysiq, anneeref, planet_type
@@ -97 +97 @@
       real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
       real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
+      real hmonsS(iip1,jjp1)
+      real summitS(iip1,jjp1)
+      real zavgS(iip1,jjp1)
       real z0S(iip1,jjp1)
 
@@ -380 +383 @@
 
       CALL datareadnc(relief,phis,alb,surfith,z0S,
-     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
+     &          zmeaS,zstdS,zsigS,zgamS,ztheS,
+     &          hmonsS,summitS,zavgS)
 
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
@@ -392 +396 @@
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
+      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,hmonsS,hmons)
+!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,summitS,summit)
+!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
 
       endif ! of if (choix_1.ne.1)
@@ -1642 +1649 @@
       call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
      .              nqtot,dtphys,real(day_ini),0.0,
-     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
+     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
+     .              hmons)
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      &              dtphys,hour_ini,