From b7f3e343f8eab7e595b45b859a055a0641afea80 Mon Sep 17 00:00:00 2001
From: Frank <franklunkeit@users.noreply.github.com>
Date: Wed, 21 Jun 2017 14:12:31 +0200
Subject: [PATCH] FL 21.06.2017:

cpl.f90: add one line diag printout (nlsg > 1)
fluxmod.f90: bug fix for entropy and energy diagnostics (missing cp in stress)
oceanmod.f90: add one line diag printout (nlsg > 1)
plasim.f90: a) bug fix for GWD entropy diagnostics (- sign)
            b) small additional print out for nprint > 0
radmod.f90: 3d co2 distribution included
---
 plasim/src/cpl.f90      |  5 +++++
 plasim/src/fluxmod.f90  |  5 +++--
 plasim/src/oceanmod.f90 |  5 +++++
 plasim/src/plasim.f90   |  6 +++++-
 plasim/src/radmod.f90   | 20 +++++++++++++-------
 5 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/plasim/src/cpl.f90 b/plasim/src/cpl.f90
index 15c25ea..5a85a81 100644
--- a/plasim/src/cpl.f90
+++ b/plasim/src/cpl.f90
@@ -486,6 +486,11 @@
 !     run lsg if its time
 !
       if (mod(kstep,naomod) == naomod-1) then
+!
+       if(nprint > 0) then
+        write(nud,*) 'in cpl: make lsg step'
+       endif
+!
        if(ncoupling == 1) then
         call cputssta(psst)
        endif
diff --git a/plasim/src/fluxmod.f90 b/plasim/src/fluxmod.f90
index 552ec47..f19b624 100644
--- a/plasim/src/fluxmod.f90
+++ b/plasim/src/fluxmod.f90
@@ -451,14 +451,15 @@
 !
       if(nentropy > 0) then
        dentropy(:,31)=zdtdt(:)/dentrot(:,NLEV)*dentrop(:)*dsigma(NLEV)  &
-     &               /ga
+     &               *acpd*(1.+adv*dentroq(:,NLEV))/ga
        if(nentro3d > 0) then
         dentro3d(:,1:NLEM,21)=0.
         dentro3d(:,NLEV,21)=dentropy(:,31)
        endif
       endif
       if(nenergy > 0) then
-       denergy(:,21)=zdtdt(:)*dp(:)/ga*dsigma(NLEV) 
+       denergy(:,21)=zdtdt(:)*acpd*(1.+adv*dq(:,NLEV))*dp(:)            &
+     &              /ga*dsigma(NLEV) 
        if(nener3d > 0) then
         dener3d(:,1:NLEM,21)=0.
         dener3d(:,NLEV,21)=denergy(:,21)
diff --git a/plasim/src/oceanmod.f90 b/plasim/src/oceanmod.f90
index c992caa..5b795f3 100644
--- a/plasim/src/oceanmod.f90
+++ b/plasim/src/oceanmod.f90
@@ -430,6 +430,11 @@
        call ntomin(nstep,ndatim(5),ndatim(4),ndatim(3),ndatim(2)        &
      &            ,ndatim(1))
        if(mypid == nroot) then
+!
+        if(nprint > 0) then
+         write(nud,*)'call for LSG coupling'
+        endif
+!         
         call clsgstep(ndatim,nstep,zsst,ztaux,ztauy,zpme,zroff,zice     &
      &               ,zheat,zfldo)
        endif
diff --git a/plasim/src/plasim.f90 b/plasim/src/plasim.f90
index 69c9883..1cfa77f 100644
--- a/plasim/src/plasim.f90
+++ b/plasim/src/plasim.f90
@@ -2788,6 +2788,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
       if(nprint > 0) then
        if(mypid==NROOT) then
           write(nud,*)'------------ SUBROUTINE GRIDPOINTD -------------'
+          write(nud,*)'NSTEP= ',nstep
+          write(nud,*)'YY,MM,DD,HH,MIN= ',ndatim(1:5)
+          write(nud,*) 'gp= ',MAXVAL(gp(1:NHOR)),MINVAL(gp(1:NHOR))
+          write(nud,*) 'dp= ',MAXVAL(dp(1:NHOR)),MINVAL(dp(1:NHOR))
        endif
        call prdbug1
       endif
@@ -2902,7 +2906,7 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
         dentropy(:,4)=dentropy(:,4)+dentro(:)
         if(nentro3d > 0) dentro3d(:,jlev,4)=dentro(:)
-        zfac=-1.*ww*tfrc(jlev)/(1.+tfrc(jlev)*delt2)
+        zfac=ww*tfrc(jlev)/(1.+tfrc(jlev)*delt2)
         dentro(:)=((du(:,jlev)+dudt(:,jlev)*deltsec2)**2                &
      &             +(dv(:,jlev)+dvdt(:,jlev)*deltsec2)**2)              &
      &             *zfac*dentrop(:)/ga*dsigma(jlev)/dentrot(:,jlev)
diff --git a/plasim/src/radmod.f90 b/plasim/src/radmod.f90
index 12dc648..791eaee 100644
--- a/plasim/src/radmod.f90
+++ b/plasim/src/radmod.f90
@@ -62,6 +62,7 @@
       real :: gmu0(NHOR)                   ! cosine of solar zenit angle
       real :: gmu1(NHOR)                   ! cosine of solar zenit angle
       real :: dqo3(NHOR,NLEV)        = 0.0 ! ozon concentration (kg/kg)
+      real :: dqco2(NHOR,NLEV)       = 0.0 ! co2 concentration (ppmv)
       real :: dtdtlwr(NHOR,NLEV)           ! lwr temperature tendencies
       real :: dtdtswr(NHOR,NLEV)           ! swr temperature tendencies
 
@@ -317,6 +318,13 @@
          dqo3cl(:,:,:) = 0.0
          call mpsurfgp('dqo3cl',dqo3cl,NHOR,NLEV*14)
       endif
+!
+!     set co2 3d-field (enable external co2 by if statement)
+!
+!
+      if(co2 > 0.) then
+       dqco2(:,:)=co2
+      endif
 !
       return
       end subroutine radini
@@ -1376,6 +1384,7 @@
       real ztau(NHOR,NLEV)      ! total transmissivity
       real zq(NHOR,NLEV)        ! modified water vapour
       real zqo3(NHOR,NLEV)      ! modified ozon
+      real zqco2(NHOR,NLEV)     ! modified co2
       real ztausf(NHOR,NLEV)    ! total transmissivity to surface
       real ztaucs(NHOR,NLEV)    ! clear sky transmissivity
       real ztaucc0(NHOR,NLEV)   ! layer transmissivity cloud
@@ -1408,8 +1417,6 @@
       zco20=0.0676*(0.01022)**0.421  ! to get a(co2)=0 for co2=0
       zh2o0a=0.846*(3.59E-5)**0.243  ! to get a(h2o)=0 for h2o=0
       zh2o0=0.832*0.0286**0.26       ! to get t(h2o)=1 for h2o=0
-      zqco2=co2*zpv2pm*1.E-6         ! co2 pp mass needed for transmissivities
-                                     ! note: co2 in ppmv not ppv
 !
 !     to make a(o3) continues at 0.01cm: 
 !
@@ -1469,7 +1476,6 @@
       zbu(:,NLEP)=zeps(:)*zst4(:,NLEP)
       zbue1(:,NLEP)=0.
       zbue2(:,NLEP)=zbu(:,NLEP)
-
 !
 !**   3) vertical loops
 !
@@ -1479,12 +1485,12 @@
       do jlev=1,NLEV
        jlep=jlev+1
        zzf1=sigma(jlev)*dsigma(jlev)/ga/100000.
-       zzf2=zfco2*zqco2*(zsigh2(jlep)-zsigh2(jlev))*0.5/ga/100000.
+       zzf2=zpv2pm*1.E-6                        !get co2 in pp mass (kg/kg-stp)
        zzf3=-1.66*acllwr*1000.*dsigma(jlev)/ga
        zsfac(:)=zzf1*zps2(:)
        zq(:,jlev)=zfh2o*zsfac(:)*dq(:,jlev)
        zqo3(:,jlev)=zfo3*zsfac(:)*dqo3(:,jlev)
-       zsumco2(:)=zzf2*zps2(:)
+       zqco2(:,jlev)=zfco2*zzf2*zsfac(:)*dqco2(:,jlev) 
        if(clgray > 0) then
         ztaucc0(:,jlev)=1.-dcc(:,jlev)*clgray
        else
@@ -1499,15 +1505,15 @@
        ztaucc(:)=1.
        zsumwv(:)=0.
        zsumo3(:)=0.
+       zsumco2(:)=0.
 !
 !     transmissivities
 !
        do jlev2=jlev,NLEV
         jlep2=jlev2+1
-        zzf2=zfco2*zqco2*(zsigh2(jlep2)-zsigh2(jlev))*0.5/ga/100000.
         zsumwv(:)=zsumwv(:)+zq(:,jlev2)
         zsumo3(:)=zsumo3(:)+zqo3(:,jlev2)
-        zsumco2(:)=zzf2*zps2(:)
+        zsumco2(:)=zsumco2(:)+zqco2(:,jlev2)
 !
 !     clear sky transmisivity
 !
-- 
GitLab