New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
tradwl.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/TRA/tradwl.F90

Last change on this file was 15651, checked in by hadjt, 2 years ago

TRA/tradwl.F90

Code to report the trends fixed. This now works without affecting the stratification (unlike the original CO6 code), but hasn't been tested in terms of the trends.

The choice of which code to report the trends too (jptra_qsr) may not be correct.
File size: 11.3 KB
Line 
1MODULE tradwl
2   !!======================================================================
3   !!                       ***  MODULE  tradwl  ***
4   !! Ocean physics: solar radiation penetration in the top ocean levels
5   !!======================================================================
6   !! History :  POLCOMS  !  1996-10  (J. Holt)  Original code
7   !!   NEMO     3.2  !  2010-03  (E. O'Dea)  Import to Nemo for use in Shelf Model
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_dwl      : trend due to the solar radiation penetration
12   !!   tra_dwl_init : solar radiation penetration initialization
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and active tracers
15   USE dom_oce         ! ocean space and time domain
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE trc_oce         ! share SMS/Ocean variables
18   USE trd_oce         ! trends: ocean variables
19   USE trdtra          ! ocean active tracers trends
20   USE in_out_manager  ! I/O manager
21   USE phycst          ! physical constants
22   USE prtctl          ! Print control
23   USE iom             ! I/O manager
24   USE fldread         ! read input fields
25   !JT
26   USE domzgr
27   USE domain
28   !JT
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_dwl        ! routine called by step.F90 (ln_tradwl=T)
33
34   !                                           !!* Namelist namtra_qsr: penetrative solar radiation
35   LOGICAL , PUBLIC ::   ln_tradwl  = .TRUE.    ! light absorption (dwl) flag
36   LOGICAL , PUBLIC ::   ln_vary_lambda  = .TRUE.    ! vary Lambda or not flag
37   
38   !! * Substitutions
39!#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
43   !! $Id$
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE tra_dwl( kt )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_qsr  ***
52      !!
53      !! ** Purpose :   Compute the temperature trend due to the solar radiation
54      !!      penetration and add it to the general temperature trend.
55      !!
56      !! ** Method  : The profile of the solar radiation within the ocean is defined
57      !!
58      !!          Jason Holt Oct 96
59      !!
60      !!          Calculates change in temperature due to penetrating
61      !!          radiation, with cooling at the surface layer
62      !!
63      !!          rad=rad0*exp(lambda*z)
64      !!
65      !!       Heat input into box is between z=K and z=K+1 is RAD(K)-RAD(K+1)
66      !!
67      !!
68      !! ** Action  : - update ta with the penetrative solar radiation trend
69      !!              - save the trend in ttrd ('key_trdtra')
70      !!
71      !!----------------------------------------------------------------------
72
73      !JT USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
74      !JT USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
75
76      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
77
78      !!
79      INTEGER, INTENT(in) ::   kt     ! ocean time-step
80      !!
81      INTEGER  ::   ji, jj, jk           ! dummy loop indices
82      INTEGER  ::   irgb                 ! temporary integers
83      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
84      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
85      !JT
86      REAL(wp), DIMENSION(jpi,jpj) ::   hbatt, qsr_tradwl
87      !JT
88      !!----------------------------------------------------------------------
89      !! HERE GO VARIABLES USED IN POLCOMS CLEAN UP LATER
90
91      integer i,j,k
92!      real*8 dtmp(n-1)
93      real*8 dtmp(jpkm1)
94      real*8 z1,z2,Rad0,Rad1,Rad2,rD,SurfOut,cp
95      logical first
96      save first
97      data first/.true./
98      !!--------------------------End of POLCOMS variables Note instead of using saves
99      !!--------------------------Could shift this into initial code
100
101      IF( kt == nit000 ) THEN
102         IF(lwp) WRITE(numout,*)
103         IF(lwp) WRITE(numout,*) 'tra_dwl : penetration of the surface solar radiation'
104         IF(lwp) WRITE(numout,*) '~~~~~~~'
105         CALL tra_dwl_init
106         IF( .NOT.ln_tradwl )   RETURN
107      ENDIF
108
109      !JT IF( l_trdtra ) THEN      ! Save ta and sa trends
110      !JT    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
111      !JT    ztrds(:,:,:) = 0.e0
112      !JT ENDIF
113
114
115      IF( l_trdtra )   THEN                  !* Save ta and sa trends
116         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
117         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
118         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
119      ENDIF
120!--------------------------------------------------------------------
121!  Set transmissivity
122!--------------------------------------------------------------------
123!
124!  Normal value
125!
126!---------------------------------------------------------------------------
127!
128!  Convert Heat fluxes to units used in POL subroutine dwl
129!
130!---------------------------------------------------------------------------
131    !cp=3986.0d0
132
133    DO jj = 2, jpj
134         DO ji = fs_2, fs_jpim1
135           qsr_tradwl(ji,jj)  = qsr(ji,jj)  * (r1_rau0_rcp)
136         ENDDO       !ji
137    ENDDO            !jj
138!--------------------------------------------------------------------------------
139 
140
141   if ( first ) then
142    do jj=2,jpjm1
143      do ji = fs_2, fs_jpim1 
144          IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
145            hbatt(ji,jj) = sum( e3t_n(ji,jj,:)*tmask(ji,jj,:) )
146        else
147            hbatt(ji,jj)= 0.
148        endif
149      enddo ! ji
150    enddo ! jj
151
152   !CALL iom_put('hbatt_tradwl', hbatt(:,:) )
153
154        rlambda2(:,:) = 0.0
155        first=.false.
156        if ( ln_vary_lambda ) then
157
158        do jj=2,jpjm1
159          do ji = fs_2, fs_jpim1   ! vector opt.
160              !IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
161
162
163              rlambda2(ji,jj)=-0.033*log(hbatt(ji,jj))+0.2583    ! JIAs formula
164              rlambda2(ji,jj)=max(0.05,rlambda2(ji,jj))     ! limit in deep water
165              rlambda2(ji,jj)=min(0.25,rlambda2(ji,jj))     ! Catch the infinities, from very shallow water/land. 10cm = 0.25
166
167            !else
168            !    rlambda2(ji,jj)= 0.25
169            !endif
170          enddo ! ji
171        enddo ! jj
172        rlambda = 0.0
173       else
174        rLambda=0.154
175       endif ! If vary lambda
176      endif ! If first
177
178      ! CALL iom_put('rlambda2_tradwl', rlambda2(:,:) )
179
180      DO jk=2,jpk
181         DO jj=2,jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183
184              IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
185
186    !--------------------------------------------------------------------
187    ! Calculate change in temperature
188    !--------------------------------------------------------------------
189    !
190    !        rad0 = hfl_in(i,j)   ! change hfl_in to qsr I assume
191
192                    rad0 = qsr_tradwl(ji,jj)
193                    rD = rLambda2(ji,jj)  +rLambda      !  Transmissivity to be used here
194                          !       if rlambda 0 then rlambda2 not zer and vica versa
195
196                    z2=gdepw_0(ji,jj,jk-1)    ! grid box is from z=z1 to z=z2
197                    z1=gdepw_0(ji,jj,jk)
198
199                    Rad2=Rad0*(exp(-z2*rD)) ! radiation entering box
200                    Rad1=Rad0*(exp(-z1*rD)) ! radiation leaving box
201
202
203                    dtmp(jk)=1.0/(e3t_0(ji,jj,jk))*(Rad2-Rad1) !change in temperature
204                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + dtmp(jk)
205                endif ! if land
206            enddo  ! ji
207         enddo  ! jj
208      enddo !jk
209
210
211      !JT IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
212      !JT    ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
213      !JT    !CEODCALL trd_mod( ztrdt, ztrds, jptra_trd_qsr, 'TRA', kt )
214      !JT ENDIF
215      !                       
216      IF( l_trdtra ) THEN      ! qsr tracers trends saved for diagnostics
217
218         !JT I think I should use jptra_qsr??
219
220         !ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
221         !ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
222         DO jk = 1, jpkm1
223            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk)
224            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk)
225         END DO
226
227         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
228         CALL trd_tra( kt, 'TRA', jp_sal, jptra_qsr, ztrds )
229         DEALLOCATE( ztrdt , ztrds )
230      ENDIF
231
232      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
233      !
234   END SUBROUTINE tra_dwl
235
236
237   SUBROUTINE tra_dwl_init
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE tra_dwl_init  ***
240      !!
241      !! ** Purpose :   Initialization for the penetrative solar radiation for Downwell routine
242      !!
243      !! ** Method  :   The profile of solar radiation within the ocean is set
244      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
245      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
246      !!      default values correspond to clear water (type I in Jerlov'
247      !!      (1968) classification.
248      !!         called by tra_qsr at the first timestep (nit000)
249      !!
250      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
251      !!
252      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
253      !!----------------------------------------------------------------------
254      INTEGER  ::   ji, jj, jk            ! dummy loop indices
255      INTEGER  ::   ios                   ! Local integer output status for namelist read
256      INTEGER  ::   irgb, ierror          ! temporary integer
257      INTEGER  ::   ioptio, nqsr          ! temporary integer
258      REAL(wp) ::   zc0  , zc1            ! temporary scalars
259      REAL(wp) ::   zc2  , zc3  , zchl    !    -         -
260      REAL(wp) ::   zsi0r, zsi1r, zcoef   !    -         -
261      !!
262      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
263      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
264      NAMELIST/namtra_dwl/  ln_tradwl, ln_vary_lambda
265      !!----------------------------------------------------------------------
266
267      REWIND( numnam_ref )            ! Read Namelist namtra_dwl in reference namelist :
268      READ  ( numnam_ref, namtra_dwl, IOSTAT = ios, ERR = 901)
269901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist')
270     
271      REWIND( numnam_cfg )            ! Read Namelist namtra_dwl in configuration namelist :
272      READ  ( numnam_cfg, namtra_dwl, IOSTAT = ios, ERR = 902)
273902   IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist')
274      !
275      IF(lwp) THEN                ! control print
276         WRITE(numout,*)
277         WRITE(numout,*) 'tra_dwl_init : '
278         WRITE(numout,*) '~~~~~~~~~~~~'
279         WRITE(numout,*) '   Namelist namtra_dwl : set the parameter of penetration'
280         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_tradwl  = ', ln_tradwl
281         WRITE(numout,*) '      Vary Lambda  (T) or not (F))             ln_vary_lambda  = ', ln_vary_lambda
282      ENDIF
283
284   END SUBROUTINE tra_dwl_init
285
286   !!======================================================================
287END MODULE tradwl
Note: See TracBrowser for help on using the repository browser.