Changeset 1208 for trunk/LMDZ.MARS/libf

Timestamp:

Mar 12, 2014, 4:30:18 PM (11 years ago)

Author:

tnavarro

Message:

added tauscaling in startfi + moved start_archive routines from dyn3d to phymars

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• aeropacity.F (modified) (4 diffs)
• ini_archive.F (moved) (moved from trunk/LMDZ.MARS/libf/dyn3d/ini_archive.F)
• lect_start_archive.F (moved) (moved from trunk/LMDZ.MARS/libf/dyn3d/lect_start_archive.F) (9 diffs)
• newstart.F (modified) (6 diffs)
• phyetat0.F90 (modified) (3 diffs)
• phyredem.F90 (modified) (3 diffs)
• physiq.F (modified) (4 diffs)
• start2archive.F (modified) (6 diffs)
• updatereffrad.F (modified) (1 diff)
• write_archive.F (moved) (moved from trunk/LMDZ.MARS/libf/dyn3d/write_archive.F)

trunk/LMDZ.MARS/libf/phymars/aeropacity.F

@@ -120 +120 @@
 c       when computing the dust opacity in each layer
 c       (this applies when doubleq and active are true)
-      INTEGER, PARAMETER :: cstdustlevel = 7
+      INTEGER, PARAMETER :: cstdustlevel0 = 7
+      INTEGER, SAVE      :: cstdustlevel
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -198 +199 @@
 !       otherwise default value read from starfi.nc file will be used)
         call getin("tauvis",tauvis)
+
+        IF (freedust) THEN
+          cstdustlevel = 1
+        ELSE
+          cstdustlevel = cstdustlevel0
+        ENDIF
+
 
         firstcall=.false.
@@ -454 +462 @@
             tauscaling(ig) = tauref(ig) * 
      &                       pplev(ig,1) / odpref / taudusttmp(ig)
-c            tauscaling(ig) = 1.e-4
         ENDDO
 
@@ -470 +477 @@
 
       IF (freedust) THEN
-       ! tauref has been initialized to 0 before.
-       DO iaer=1,naerdust
-        DO l=1,nlayer
-         DO ig=1,ngrid
-          tauref(ig) = tauref(ig) + 
-     &                aerosol(ig,l,iaerdust(iaer))
-          ENDDO
-        ENDDO
-       ENDDO
+        ! tauref has been initialized to 0 before.
+        DO iaer=1,naerdust
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              tauref(ig) = tauref(ig) + 
+     &                    aerosol(ig,l,iaerdust(iaer))
+            ENDDO
+          ENDDO
+        ENDDO
+        tauref(:) = tauref(:) * odpref / pplev(:,1)
       ENDIF
       

trunk/LMDZ.MARS/libf/phymars/lect_start_archive.F

@@ -2 +2 @@
      &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,co2ice,h,phisold_newgrid,
-     &     q,qsurf,surfith,nid)
+     &     q,qsurf,tauscaling,surfith,nid)
 c=======================================================================
 c
@@ -109 +109 @@
       REAL emis(ngrid)
       REAL q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
+      REAL tauscaling(ngrid) ! dust conversion factor
 c     REAL phisfi(ngrid)
 
@@ -136 +137 @@
       real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
+      real tauscalingS(iip1,jjp1)
 
       real ptotal, co2icetotal
@@ -163 +165 @@
       real, dimension(:,:), allocatable :: emisold
       real, dimension(:,:,:,:), allocatable :: qold
+      real, dimension(:,:), allocatable :: tauscalingold
 
       real tab_cntrl(100)
@@ -337 +340 @@
       allocate(mlayerold(nsoilold))
       allocate(qsurfold(imold+1,jmold+1,nqtot))
+      allocate(tauscalingold(imold+1,jmold+1))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -662 +666 @@
 
 c-----------------------------------------------------------------------
-c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
+c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf, tauscaling)
 c-----------------------------------------------------------------------
  
@@ -742 +746 @@
          PRINT*, "lect_start_archive: Failed loading <q2surf>"
          CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "tauscaling", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: <tauscaling> not in file"
+         tauscalingold(:,:) = -1
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tauscalingold)
+#else
+        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tauscalingold)
+#endif
+        IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "lect_start_archive: Failed loading <tauscaling>"
+           PRINT*, NF_STRERROR(ierr)
+           CALL abort
+        ENDIF
       ENDIF
 c
@@ -1019 +1040 @@
       call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
+c Dust conversion factor
+      call interp_horiz (tauscalingold,tauscalings,imold,jmold,iim,jjm,
+     &                   1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tauscalings,tauscaling)
 c     write(46,*) 'emis',emis
 c-----------------------------------------------------------------------
@@ -1378 +1402 @@
       deallocate(mlayerold)
       deallocate(qsurfold)
+      deallocate(tauscalingold)
       deallocate(var,varp1)
 

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

@@ -109 +109 @@
       REAL co2ice(ngridmx)      ! CO2 ice layer
       REAL emis(ngridmx)        ! surface emissivity
+      REAL tauscaling(ngridmx) ! dust conversion factor
+      REAL tauscadyn(iip1,jjp1) ! dust conversion factor on the dynamics grid
       REAL,ALLOCATABLE :: qsurf(:,:)
       REAL q2(ngridmx,nlayermx+1)
@@ -426 +428 @@
      &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,
-     &   surfith,nid)
+     &   tauscaling,surfith,nid)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -447 +449 @@
         CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      .        day_ini,time,
-     .        tsurf,tsoil,emis,q2,qsurf,co2ice)
+     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
         
         ! copy albedo and soil thermal inertia
@@ -858 +860 @@
              endif
              
-c       q=profile : initialize tracer with a given profile
+c       convert dust from virtual to true values
 c       --------------------------------------------------
         else if (trim(modif) .eq. 'freedust') then
+         if (minval(tauscaling) .lt. 0) then
+           write(*,*) 'WARNING conversion factor negative'
+           write(*,*) 'This is probably because it was not present
+     &in the file'
+           write(*,*) 'A constant conversion is used instead.'
+           tauscaling(:) = 1.e-3
+         endif
+         CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscadyn)
           do l=1,llm
             do j=1,jjp1
               do i=1,iip1
                 if (igcm_dust_number .ne. 0) 
-     &          q(i,j,l,igcm_dust_number)=
-     &                q(i,j,l,igcm_dust_number) * 1e-3 ! grosso modo
+     &            q(i,j,l,igcm_dust_number) =
+     &            q(i,j,l,igcm_dust_number) * tauscadyn(i,j)
                 if (igcm_dust_mass .ne. 0) 
-     &          q(i,j,l,igcm_dust_mass)=
-     &                q(i,j,l,igcm_dust_mass)   * 1e-3 ! grosso modo
+     &            q(i,j,l,igcm_dust_mass) =
+     &            q(i,j,l,igcm_dust_mass) * tauscadyn(i,j)
                 if (igcm_ccn_number .ne. 0) 
-     &          q(i,j,l,igcm_ccn_number)=
-     &                q(i,j,l,igcm_ccn_number) * 1e-3 ! grosso modo
+     &            q(i,j,l,igcm_ccn_number) =
+     &            q(i,j,l,igcm_ccn_number) * tauscadyn(i,j)
                 if (igcm_ccn_mass .ne. 0) 
-     &          q(i,j,l,igcm_ccn_mass)=
-     &                q(i,j,l,igcm_ccn_mass)   * 1e-3 ! grosso modo
+     &            q(i,j,l,igcm_ccn_mass) =
+     &            q(i,j,l,igcm_ccn_mass) * tauscadyn(i,j)
               end do
             end do
           end do
+
+          tauscaling(:) = 1.
 
          ! We want to have the very same value at lon -180 and lon 180
@@ -888 +900 @@
              end do
           end do
+
+          write(*,*) 'done rescaling to true vale'
 
 c       ini_q : Initialize tracers for chemistry
@@ -1529 +1543 @@
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      .              dtphys,hour_ini,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf)
+     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
 
 c=======================================================================

trunk/LMDZ.MARS/libf/phymars/phyetat0.F90

@@ -1 +1 @@
 subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
-                     day_ini,time0,tsurf,tsoil,emis,q2,qsurf,co2ice)
-
+                     day_ini,time0,tsurf,tsoil,emis,q2,qsurf,co2ice, &
+                     tauscaling)
 !  use netcdf
   use infotrac, only: nqtot, tname
@@ -57 +57 @@
   real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
   real,intent(out) :: co2ice(ngrid) ! co2 ice cover
+  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
 
 !======================================================================
@@ -236 +237 @@
 
 
+! Dust conversion factor
+call get_field("tauscaling",tauscaling,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <tauscaling> not in file"
+  tauscaling(:) = -1
+else
+  write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
+             minval(tauscaling), maxval(tauscaling)
+endif
+
+
 ! Surface temperature :
 call get_field("tsurf",tsurf,found,indextime)

trunk/LMDZ.MARS/libf/phymars/phyredem.F90

@@ -143 +143 @@
 
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
-                    phystep,time,tsurf,tsoil,co2ice,emis,q2,qsurf)
+                    phystep,time,tsurf,tsoil,co2ice,emis,q2,qsurf,&
+                    tauscaling)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -162 +163 @@
   real,intent(in) :: q2(ngrid,nlay+1)
   real,intent(in) :: qsurf(ngrid,nq)
+  real,intent(in) :: tauscaling(ngrid)
   
   integer :: iq
@@ -192 +194 @@
   call put_field("q2","pbl wind variance",q2,time)
   
+  ! Dust conversion factor
+  ! Only to be read by newstart to convert to actual dust values
+  ! Or by any user who wants to reconstruct dust, opacity from the start files.
+  call put_field("tauscaling","dust conversion factor",tauscaling,time)
+
   ! Tracers on the surface
   ! preliminary stuff: look for water vapour & water ice tracers (if any)

trunk/LMDZ.MARS/libf/phymars/physiq.F

@@ -436 +436 @@
      &         nsoilmx,ngrid,nlayer,nq,
      &         day_ini,time_phys,
-     &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
 #else
 #include "meso_inc/meso_inc_ini.F"
@@ -1582 +1582 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf,tsoil,co2ice,emis,q2,qsurf)
+     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
           
          ENDIF
@@ -2136 +2136 @@
 c            call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
 c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
+c             call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
+c     &                       'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
+c             call WRITEDIAGFI(ngrid,'subice','sublimated ice',
+c     &                       'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
              call WRITEDIAGFI(ngrid,'dqsdust',
      &                        'deposited surface dust mass',
@@ -2148 +2152 @@
              call WRITEDIAGFI(ngrid,'dustN','Dust number',
      &                        'part/kg',3,ndust)
+c             call WRITEDIAGFI(ngrid,"tauscaling",
+c     &                    "dust conversion factor"," ",2,tauscaling)
 #ifdef MESOINI
 !     !!! to initialize mesoscale we need scaled variables

trunk/LMDZ.MARS/libf/phymars/start2archive.F

@@ -67 +67 @@
       REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
       REAL co2ice(ngridmx)      ! CO2 ice layer
+      REAL tauscaling(ngridmx) ! dust conversion factor
       REAL q2(ngridmx,llm+1)
       REAL,ALLOCATABLE :: qsurf(:,:)
@@ -82 +83 @@
       REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
       REAL co2iceS(ip1jmp1)
+      REAL tauscalingS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
       REAL,ALLOCATABLE :: qsurfS(:,:)
@@ -147 +149 @@
 
       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
-     &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice)
+     &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice,
+     &      tauscaling)
 
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
@@ -242 +245 @@
 c q2 --> q2S
 c qsurf --> qsurfS
+c tauscaling --> tauscalingS
 c
 c-----------------------------------------------------------------------
@@ -253 +257 @@
       call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
       call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscalingS)
 
 c=======================================================================
@@ -341 +346 @@
       call write_archive(nid,ntime,'co2ice','couche de glace co2',
      &  'kg/m2',2,co2iceS)
+      call write_archive(nid,ntime,'tauscaling',
+     &  'dust conversion factor',' ',2,tauscalingS)
       call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
       call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)

trunk/LMDZ.MARS/libf/phymars/updatereffrad.F

@@ -128 +128 @@
         IF (water.AND.activice) THEN 
          IF (microphys) THEN
-         
-          IF (firstcall) THEN
-            DO l=1,nlayer
-              DO ig=1,ngrid
-                call updaterice_micro(pq(ig,l,igcm_h2o_ice),
-     &                                pq(ig,l,igcm_ccn_mass),
-     &                                pq(ig,l,igcm_ccn_number),
-     &                                1.e-3,rice(ig,l),
-     &                                rhocloud(ig,l))
-                nuice(ig,l) = nuice_ref
-              ENDDO
-            ENDDO
-          firstcall = .false.
+        
 c    At firstcall, the true number and true mass of cloud condensation nuclei are not known.
 c    Indeed it is scaled on the prescribed dust opacity via a 'tauscaling' coefficient
-c    computed after radiative transfer.
-          ELSE
-            DO l=1,nlayer
-              DO ig=1,ngrid
-                call updaterice_micro(pq(ig,l,igcm_h2o_ice),
-     &                                pq(ig,l,igcm_ccn_mass),
-     &                                pq(ig,l,igcm_ccn_number),
-     &                                tauscaling(ig),rice(ig,l),
-     &                                rhocloud(ig,l))
-                nuice(ig,l) = nuice_ref
-               ENDDO
-             ENDDO
-          ENDIF ! of if firstcall
+c    computed after radiative transfer. If tauscaling is not in startfi, we make an assumption for its value.
+
+          IF (firstcall) THEN
+            IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
+            IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
+            firstcall = .false.
+          ENDIF
+ 
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              call updaterice_micro(pq(ig,l,igcm_h2o_ice),
+     &                              pq(ig,l,igcm_ccn_mass),
+     &                              pq(ig,l,igcm_ccn_number),
+     &                              tauscaling(ig),rice(ig,l),
+     &                              rhocloud(ig,l))
+              nuice(ig,l) = nuice_ref
+            ENDDO
+          ENDDO
           
         ELSE ! if not microphys