Changeset 1375

Timestamp:

Feb 11, 2015, 10:31:46 PM (10 years ago)

Author:

aslmd

Message:

modification to physiq for dust bomb simulation

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• aeropacity.F (modified) (5 diffs)
• physiq.F (modified) (4 diffs)

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

@@ -9 +9 @@
      &                      nqdust
       use comgeomfi_h, only: lati, sinlat ! grid point latitudes (rad)
+#ifdef DUSTSTORM
+      use comgeomfi_h, only: long
+      use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
+#endif
       use planete_h
       USE comcstfi_h
@@ -108 +112 @@
       REAL topdust0(ngrid)
 
+#ifdef DUSTSTORM
+!! Local dust storms
+      logical localstorm        ! =true to create a local dust storm
+      real taulocref,ztoploc,radloc,lonloc,latloc  ! local dust storm parameters
+      real reffstorm, yeah
+      REAL ray(ngrid) ! distance from dust storm center
+      REAL tauuser(ngrid) ! opacity perturbation due to dust storm
+      REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
+      REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
+      real l_top ! layer of the storm's top
+      REAL zalt(ngrid, nlayer) ! useful factor to compute l_top 
+      REAL zdqnorm(ngrid, nlayer,2)
+#endif
+
 c   local saved variables
 c   ---------------------
@@ -189 +207 @@
 
 
+#ifndef DUSTSTORM
         firstcall=.false.
+#endif
 
       END IF
@@ -437 +457 @@
 c -----------------------------------------------------------------
 
+#ifdef DUSTSTORM
+c -----------------------------------------------------------------
+c the quantity of dust to add at the first time step is calculated
+c to match a tunable opacity perturbation.
+c -----------------------------------------------------------------
+      IF (firstcall) THEN
+        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)
+c--------------------------------------------------
+c  Parameters of the opacity perturbation
+c--------------------------------------------------
+
+      iaer=1  !!!! PROVISOIRE !!!!
+
+        write(*,*) "Add a local storm ?"
+        localstorm=.true. ! default value
+        call getin("localstorm",localstorm)
+        write(*,*) " localstorm = ",localstorm
+
+        IF (localstorm) THEN
+          WRITE(*,*) "********************"
+          WRITE(*,*) "ADDING A LOCAL STORM"
+          WRITE(*,*) "********************"
+
+          write(*,*) "ref opacity of local dust storm"
+              taulocref = 4.25 ! default value
+              call getin("taulocref",taulocref)
+              write(*,*) " taulocref = ",taulocref
+
+          write(*,*) "target altitude of local storm (km)"
+              ztoploc = 10.0 ! default value
+              call getin("ztoploc",ztoploc)
+              write(*,*) " ztoploc = ",ztoploc
+
+          write(*,*) "radius of dust storm (degree)"
+              radloc = 0.5 ! default value
+              call getin("radloc",radloc)
+              write(*,*) " radloc = ",radloc
+
+          write(*,*) "center longitude of storm (deg)"
+              lonloc = 25.0 ! default value
+              call getin("lonloc",lonloc)
+              write(*,*) " lonloc = ",lonloc
+
+          write(*,*) "center latitude of storm (deg)"
+              latloc = -2.5 ! default value
+              call getin("latloc",latloc)
+              write(*,*) " latloc = ",latloc
+        
+          write(*,*) "reff storm (mic) 0. for background"
+              reffstorm = 0.0 ! default value
+              call getin("reffstorm",reffstorm)
+              write(*,*) " reffstorm = ",reffstorm
+
+      DO ig=1,ngrid
+
+
+c---------------------------------------
+c        distance to the center:
+c-----------------------------------------
+
+      ray(ig)=SQRT((lati(ig)*180./pi-latloc)**2 +
+     &          (long(ig)*180./pi -lonloc)**2)
+
+
+      !! transition factor for storm
+      !! increase factor ray diff for steepness
+      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
+
+c-------------------------------------------------
+c           Tau's new map:
+c------------------------------------------
+
+      tauuser(ig)=max(tauref(ig) * pplev(ig,1) /odpref  , 
+     &          taulocref * yeah)
+
+c---------------------------------------------------------
+c           compute l_top
+c----------------------------------------------------------
+
+          DO l=1,nlayer
+            zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
+     &                      / g / 44.01
+     &                    * 8.31 * 210.
+                IF (     (ztoploc .lt. zalt(ig,l)  )
+     &          .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
+          ENDDO
+
+c---------------------------------------------------------
+c           change reffrad if ever needed
+c----------------------------------------------------------
+
+      IF (reffstorm .gt. 0.) THEN
+          DO l=1,nlayer
+             IF (l .lt. l_top+1) THEN
+                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
+     &          * 1.e-6 * yeah )
+             ENDIF
+          ENDDO
+      ENDIF
+
+
+c---------------------------------------------------------------------------
+      ENDDO !! boucle sur ngridmx
+c---------------------------------------------------------------------------
+c----------------------------------------------------------------------------
+c    compute int_factor
+c---------------------------------------------------------------------------
+
+      DO ig=1,ngrid
+          int_factor(ig)=0.
+          DO l=1,nlayer
+             IF (l .lt. l_top+1) THEN
+                      int_factor(ig) =
+     &                int_factor(ig) +
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+             ENDIF
+          ENDDO
+          DO l=1, nlayer
+
+c-------------------------------------------------------------------------
+c Mass mixing ratio perturbation due to the local dust storm in each
+c layer
+c-------------------------------------------------------------------------
+          more_dust(ig,l,1)=
+     &                     (tauuser(ig)-(tauref(ig)
+     &                      * pplev(ig,1) /odpref)) /
+     &                      int_factor(ig)
+          more_dust(ig,l,2)=
+     &                     (tauuser(ig)-(tauref(ig) *
+     &                     pplev(ig,1) /odpref))
+     &                      / int_factor(ig) *
+     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
+     &                      * r3n_q 
+          ENDDO
+      ENDDO
+
+c--------------------------------------------------------------------------------------
+c   quantity of dust which is added at the first time step 
+c     in dynamical core.
+c--------------------------------------------------------------------------------------
+
+        DO l=1, l_top
+          zdqnorm(:,l,1) = more_dust(:,l,1)
+          zdqnorm(:,l,2) = more_dust(:,l,2)
+          pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)+ zdqnorm(:,l,1)
+          pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)+ zdqnorm(:,l,2)
+        ENDDO
+        ENDIF
+        tauref(:)=0.
+      ENDIF  !firstcall
+#endif
+
       IF (freedust) THEN
           tauscaling(:) = 1.
@@ -512 +694 @@
         end do
       ENDDO
+
+#ifdef DUSTSTORM
+      IF (firstcall) THEN
+        firstcall=.false.
+      ENDIF
+#endif
+
 c -----------------------------------------------------------------
 c Density scaled opacity and column opacity output

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

@@ -192 +192 @@
       INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
       INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
-     
+#ifdef DUSTSTORM
+      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
+#endif
 c     Variables used by the water ice microphysical scheme:
       REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
@@ -516 +518 @@
       zdtsurf(:)=0
       dqsurf(:,:)=0
+      pq_tmp(:,:,:)=0
       igout=ngrid/2+1 
 
@@ -667 +670 @@
 c          Radiative transfer
 c          ------------------
+
+#ifdef DUSTSTORM
+           pq_tmp(1:ngrid,1:nlayer,1)=
+     & pq(1:ngrid,1:nlayer,igcm_dust_mass)
+
+           pq_tmp(1:ngrid,1:nlayer,2)=
+     & pq(1:ngrid,1:nlayer,igcm_dust_number)
+         
+#endif
+
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
      $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
@@ -672 +685 @@
      $     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
      $     nuice,co2ice)
+
+#ifdef DUSTSTORM
+           pdq(1:ngrid,1:nlayer,igcm_dust_mass)= 
+     & pdq(1:ngrid,1:nlayer,igcm_dust_mass)
+     & - pq_tmp(1:ngrid,1:nlayer,1)
+     & +pq(1:ngrid,1:nlayer,igcm_dust_mass)
+
+
+           pdq(1:ngrid,1:nlayer,igcm_dust_number)=
+     & pdq(1:ngrid,1:nlayer,igcm_dust_number)
+     & - pq_tmp(1:ngrid,1:nlayer,2)
+     & +pq(1:ngrid,1:nlayer,igcm_dust_number)
+#endif
 
 c          Outputs for basic check (middle of domain)