Changeset 486 for trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

Timestamp:

Dec 21, 2011, 10:09:07 AM (13 years ago)

Author:

rwordsworth

Message:

Variable ice timestep added for ice evolution algorithm.
Sedimentation improved to allow consistent dust (L. Kerber).
Bug linked to allocatable matrices removed from rcm1d.F.
Treatment of initial aerosol radii in callcorrk.F90 improved.

File:

• trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 (modified) (16 diffs)

trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

@@ -211 +211 @@
 !     Effective radius and variance of the aerosols
 
-!     CO2 ice:
-         DO l=1,nlayer
-            DO ig=1,ngrid
-               reffrad(ig,l,1)  = 1.e-4
-               nueffrad(ig,l,1) = 0.1 
+         do iaer=1,naerkind
 !     these values will change once the microphysics gets to work
 !     UNLESS tracer=.false., in which case we should be working with
 !     a fixed aerosol layer, and be able to define reffrad in a 
 !     .def file. To be improved!
-            ENDDO
-         ENDDO
-
-!     H2O ice:
-         if(naerkind.eq.2)then
-            DO l=1,nlayer
-               DO ig=1,ngrid
-                  reffrad(ig,l,naerkind)  = 1.e-5
-                  nueffrad(ig,l,naerkind) = 0.1
-               ENDDO  
-            ENDDO
-         endif
+
+            if(iaer.eq.1)then ! CO2 ice
+               do l=1,nlayer
+                  do ig=1,ngrid
+                     reffrad(ig,l,iaer)  = 1.e-4
+                     nueffrad(ig,l,iaer) = 0.1 
+                  enddo
+               enddo
+            endif
+
+            if(iaer.eq.2)then ! H2O ice
+               do l=1,nlayer
+                  do ig=1,ngrid
+                     reffrad(ig,l,iaer)  = 1.e-5
+                     nueffrad(ig,l,iaer) = 0.1 
+                  enddo
+               enddo
+            endif
+
+            if(iaer.eq.3)then ! dust
+               do l=1,nlayer
+                  do ig=1,ngrid
+                     reffrad(ig,l,iaer)  = 1.e-5
+                     nueffrad(ig,l,iaer) = 0.1 
+                  enddo
+               enddo
+            endif
+ 
+            if(iaer.gt.3)then
+               print*,'Error in callcorrk, naerkind is too high.'
+               print*,'The code still needs generalisation to arbitrary'
+               print*,'aerosol kinds and number.'
+               call abort
+            endif
+
+         enddo
 
          print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
@@ -248 +268 @@
          Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
 
-
          if((igcm_h2o_vap.eq.0) .and. varactive)then
             print*,'varactive in callcorrk but no h2o_vap tracer.'
@@ -268 +287 @@
          enddo
       enddo
-
 
       if(kastprof)then
@@ -426 +444 @@
                end do
 
-
 !     longwave
                DO nw=1,L_NSPECTI 
@@ -465 +482 @@
                end do
             end do
-
 
             ! test / correct for freaky s. s. albedo values
@@ -515 +531 @@
             !tauaero(L_LEVELS+1,iaer) = 0.
          end do
-         !print*,'Note changed tauaero BCs in callcorrk!'
 
 !     Albedo and emissivity
@@ -533 +548 @@
          acosz=mu0(ig)          ! cosine of sun incident angle
       endif
-
 
 !-----------------------------------------------------------------------
@@ -591 +605 @@
       end do
 
-
-
 !-----------------------------------------------------------------------
 !     kcm mode only
       if(kastprof)then
 
+         ! initial values equivalent to mugaz
          DO l=1,nlayer
             muvarrad(2*l)   = mugaz
@@ -602 +615 @@
          END DO
 
-         do k=1,L_LEVELS
-            qvar(k) = 0.0
-         end do
-         print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
+         !do k=1,L_LEVELS
+         !   qvar(k) = 0.0
+         !end do
+         !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
       endif
 
@@ -633 +646 @@
          qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
 
-
-      endif
-
+      endif
 
       ! Keep values inside limits for which we have radiative transfer coefficients
@@ -678 +689 @@
       tmid(L_LEVELS) = tlevrad(L_LEVELS)
 
-
       ! test for out-of-bounds pressure
       if(plevrad(3).lt.pgasmin)then
@@ -695 +705 @@
             print*,'Minimum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
-            print*,'WARNING, OVERRIDING FOR TEST'
-            !call abort
+            !print*,'WARNING, OVERRIDING FOR TEST'
+            call abort
          elseif(tlevrad(k).gt.tgasmax)then
             print*,'Maximum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
-            print*,'WARNING, OVERRIDING FOR TEST'
-            !print*, 'T=',pt
-            !call abort
+            !print*,'WARNING, OVERRIDING FOR TEST'
+            call abort
          endif
       enddo
@@ -735 +744 @@
                  tmid,pmid,taugsurf,qvar,muvarrad)
 
-
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
                  acosz,stel_fract,gweight,                         &
                  nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
-                 !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
 
@@ -753 +760 @@
          end if
 
-
-
-
 !-----------------------------------------------------------------------
 !     Longwave
@@ -846 +850 @@
       if(specOLR)then
         if(ngrid.ne.1)then
-          !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
-          !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
+          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
+          call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
         endif
       endif
@@ -892 +896 @@
       endif
 
-      !!! see physiq.F for explanations about CLFvarying. This is temporary.
+      ! see physiq.F for explanations about CLFvarying. This is temporary.
       if (lastcall .and. .not.CLFvarying) then
         IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )