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.
flx_core.h90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/flx_core.h90 @ 812

Last change on this file since 812 was 812, checked in by ctlod, 16 years ago

missing continuation sign & in the computation of dqla_ice(:,:), see ticket #66

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.9 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                    ***  flx_core.h90  ***
3   !!----------------------------------------------------------------------
4   !!   flx     : update surface thermohaline fluxes from the NCAR data set
5   !!             read in a NetCDF file
6   !!----------------------------------------------------------------------
7   !! * Modules used     C A U T I O N  already defined in flxmod.F90
8
9   !! * Shared module variables
10   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
11      gsss,            &  !: SSS mean on nfbulk ocean time step
12      gu  , gv ,       &  !: (un,vn)(jk=1) mean over nfbulk ocean time step
13      !                   ! these variables are used for output in diawri
14      qlw_oce  ,       &  !: long wave flx over ocean
15      qla_oce  ,       &  !: latent heat flx over ocean
16      qsb_oce  ,       &  !: sensible heat flx over ocean
17      qlw_ice  ,       &  !: long wave flx over ice
18      qsb_ice             !: sensible heat flx over ice
19
20   !! * Module variables
21   INTEGER, PARAMETER ::   jpfile  = 8  ! maximum number of files to read
22   CHARACTER (LEN=32), DIMENSION (jpfile) ::  clvarname 
23   CHARACTER (LEN=50), DIMENSION (jpfile) ::  clname
24
25   INTEGER   ::      isnap
26   INTEGER, DIMENSION(jpfile) ::   &
27      numflxall,           &  ! logical units for surface fluxes data
28      nrecflx  , nrecflx2     ! first and second record to be read in flux file
29   REAL(wp), DIMENSION(jpfile) ::   &
30      freqh                 ! Frequency for each forcing file (hours)
31      !                     ! a negative value of -12 corresponds to monthly
32   REAL(wp), DIMENSION(jpi,jpj,jpfile,2) ::   &
33      flxdta                !:  set of NCAR 6hourly/daily/monthly fluxes
34   LOGICAL  :: &
35        &   ln_2m   = .FALSE.     !: logical flag for height of air temp. and hum.
36   REAL(wp) :: alpha_precip=1.    !: multiplication factor for precipitation
37   !!----------------------------------------------------------------------
38   !! History :
39   !!   9.0  !  04-08  (U. Schweckendiek)  Original code
40   !!        !  05-04  (L. Brodeau, A.M. Treguier)  severals modifications:
41   !!        !         - new bulk routine for efficiency
42   !!        !         - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
43   !!        !         - file names and file characteristics in namelist
44   !!        !         - Implement reading of 6-hourly fields   
45   !!        !  06-02  (S. Masson, G. Madec)  IOM read + NEMO "style"
46   !!        !  12-06  (L. Brodeau) handle both 2m and 10m input fields
47   !!----------------------------------------------------------------------
48   !!   OPA 9.0 , LOCEAN-IPSL (2006)
49   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/SBC/flx_core.h90,v 1.5 2007/06/05 10:41:37 opalod Exp $
50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE flx( kt )
56      !!---------------------------------------------------------------------
57      !!                     ***  ROUTINE flx  ***
58      !!                   
59      !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
60      !!                and momentum fluxes (tau)
61      !!      to the ocean at each time step.
62      !!
63      !! ** Method  :   Read NCAR data in a NetCDF file
64      !!               (file names and frequency of inputs specified in namelist)
65      !!          precipitation:  total (rain+snow)
66      !!          precipitation:  snow only
67      !!          u10,v10 -> scalar wind at 10m in m/s -  ON 'T' GRID POINTS!!!
68      !!          solar radiation (short wave)  in W/m2
69      !!          thermal radiation (long wave) in W/m2
70      !!          specific humidity             in %
71      !!          temperature at 10m            in degrees K
72      !!
73      !!  ** Action  :   
74      !!          call flx_blk_albedo to compute ocean and ice albedo
75      !!          Calculates forcing fluxes to input into ice model
76      !          or to be used directly for the ocean in ocesbc.
77      !!
78      !!  ** Outputs
79      !!          COMMENTS TO BE ADDED -units to be verified!
80      !!          taux    : zonal      wind stress on "u" points  (N/m2)
81      !!          tauy    : meridional wind stress on "v" points  (N/m2)
82      !!          qsr_oce : Solar flux over the ocean (W/m2)
83      !!          qnsr_oce: longwave flux over the ocean (W/m2)
84      !!          qsr_ice : solar flux over the ice (W/m2)
85      !!          qnsr_ice: longwave flux over the ice (W/m2)
86      !!          qla_ice : latent flux over sea-ice (W/m2)
87      !!          dqns_ice: total heat fluxes sensitivity over ice (dQ/dT)
88      !!          dqla_ice: latent heat flux sensitivity over ice (dQla/dT)
89      !!          fr1_i0
90      !!          fr2_i0
91      !!          tprecip
92      !!          sprecip
93      !!          evap
94      !!         
95      !! caution : now, in the opa global model, the net upward water flux is
96      !! -------   with mm/day unit,i.e. Kg/m2/s.
97      !!---------------------------------------------------------------------
98      !! * modules used
99      USE iom
100      USE restart         
101      USE par_oce
102      USE flx_oce
103      USE blk_oce         ! bulk variable
104      USE taumod
105      USE bulk
106      USE phycst
107      USE lbclnk
108      USE dtatem          ! FOR Ce = F(SST(levitus)):
109
110      !! * arguments
111      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
112
113      !! * Local declarations
114      INTEGER , PARAMETER :: jpmaxtime = 366*4   ! maximum time steps in file will be for a 6 hourly field
115                                                 ! and a leap year (necessary variable for flinopen??)
116      INTEGER ::  irectot, irecflx
117      INTEGER ::  ihour                ! current hour of the day at which we read the forcings
118      INTEGER ::   &
119         imois, imois2,             &  ! temporary integers
120         i15  , iman  , inum ,      &  !    "          "
121         ifpr , jfpr  ,             &  ! frequency of index for print in prihre
122         jj   , ji                     !  Loop indices
123      REAL(wp) ::   &
124         zadatrjb,                  &  ! date (noninteger day) at which we read the forcings
125         zxy    ,                   &  ! scalar for temporal interpolation
126         zcof , zzu , zzv              ! scalar
127      REAL(wp), DIMENSION(jpi,jpj, jpfile) ::   &
128         flxnow             ! input flux values at current time step
129      REAL(wp), DIMENSION(jpi,jpj) ::   &
130         dqlw_ice ,      &  ! long wave heat flx sensitivity over ice
131         dqsb_ice ,      &  ! sensible heat flx sensitivity over ice
132         alb_oce_os,     &  ! albedo of the ocean under overcast sky
133         alb_ice_os,     &  ! albedo of the ice under overcast sky
134         alb_ice_cs,     &  ! albedo of ice under clear sky
135         alb_oce_cs,     &  ! albedo of the ocean under clear sky
136         zsst,           &  ! nfbulk : mean SST
137         zsss,           &  ! nfbulk : mean tn_ice(:,:,1)
138         zut,            &  ! nfbulk : mean U at T-point
139         zvt,            &  ! nfbulk : mean V at T-point
140         dUnormt,        &  ! scalar wind (norm) on T points
141         tauxt, tauyt,   &  ! wind stress computed at T-point
142         qsatw,          &  ! specific humidity at zSST
143         qsat,           &  ! specific humidity at zsss
144         Ch,             &  ! coefficient for sensible heat flux
145         Ce,             &  ! coefficient for latent heat flux
146         Cd,             &  ! coefficient for wind stress
147         zt_zu,          &  !: air temperature at wind speed height
148         zq_zu              !: air spec. hum.  at wind speed height
149      REAL(wp), PARAMETER ::   &
150         rhoa  = 1.22,    &  ! air density
151         cp    = 1000.5,  &  ! specific heat of air
152         Lv    = 2.5e6,   &  ! latent heat of vaporization
153         Ls    = 2.839e6, &  ! latent heat of sublimation
154         Stef  = 5.67e-8, &  ! Stefan Boltzmann constant
155         Cice = 1.63e-3      ! transfert coefficient over ice
156      REAL(wp), DIMENSION(jpi,jpj) ::   &
157         catm1               !
158
159      NAMELIST /namcore/  clname, clvarname, freqh, ln_2m, alpha_precip
160      !!---------------------------------------------------------------------
161
162      !! =============================================== !!
163      !!   PART A:   READING FLUX FILES WHEN NECESSARY   !!
164      !! =============================================== !!
165
166      !--- calculation for monthly data
167      ! i15 = 0 if first half of current month
168      ! i15 = 1 if second half of current month
169      i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
170      iman  = 12
171      imois = nmonth + i15 - 1
172      IF( imois == 0 ) imois = iman
173      imois2 = nmonth
174
175      !--- calculation for 6-hourly data
176      !    we need the fraction of day, measured in hours at the date of forcings.
177      !    This is the time step before the date we are calculating, zadatrj
178      !    "adatrj" (real time in days (noninteger))
179      zadatrjb = adatrj - rdt / rday
180      ihour = INT( ( zadatrjb - INT(zadatrjb) ) * 24 )
181
182      !                                         ! ------------------------ !
183      IF( kt == nit000 ) THEN                   !   first call kt=nit000   !
184         !                                      ! ------------------------ !
185         CALL core_rst ( kt, 'READ' )
186
187         !--- Initializes default values of file names, frequency of forcings
188         !    and variable to be read in the files, before reading namelist
189         clname(1) = 'precip.nc'   ; freqh(1) = -12 ; clvarname(1) = 'precip'     ! monthly
190         clname(2) = 'u10.nc'      ; freqh(2) =   6 ; clvarname(2) = 'u10'        ! 6 hourly         
191         clname(3) = 'v10.nc'      ; freqh(3) =   6 ; clvarname(3) = 'v10'        ! 6 hourly   
192         clname(4) = 'q10.nc'      ; freqh(4) =   6 ; clvarname(4) = 'q10'        ! 6 hourly   
193         clname(5) = 'radsw.nc'    ; freqh(5) =  24 ; clvarname(5) = 'radsw'      ! daily
194         clname(6) = 'radlw.nc'    ; freqh(6) =  24 ; clvarname(6) = 'radlw'      ! daily
195         clname(7) = 't10.nc'      ; freqh(7) =   6 ; clvarname(7) = 't10'        ! 6 hourly
196         clname(8) = 'snow.nc'     ; freqh(8) = -12 ; clvarname(8) = 'snow'       ! monthly
197
198         !--- Read Namelist namcore : OMIP/CORE
199         REWIND ( numnam )
200         READ   ( numnam, namcore )
201
202         ! in case of ln_2m, air temp. and humidity are given at 2 m, thus file name and variable name are changed
203         IF ( ln_2m ) THEN
204           clname(4) = 'q2.nc'      ; clvarname(4) = 'q2'     
205           clname(7) = 't2.nc'      ; clvarname(7) = 't2'   
206           ! but if for some reason, clname and varname were also in the name list, we screw them up !
207           ! re-read the namelist, to restore clname, varname to their namelist value
208           REWIND ( numnam )
209           READ   ( numnam, namcore )
210         ENDIF
211
212         IF(lwp) THEN
213            WRITE(numout,*)' '
214            WRITE(numout,*)' flx_core : global CORE fields in NetCDF format'
215            WRITE(numout,*)' ~~~~~~~~~~ '
216            WRITE(numout,*) '           ln_2m        = ', ln_2m
217            WRITE(numout,*) '           alpha_precip = ', alpha_precip
218            WRITE(numout,*) '           list of files and frequency (hour), or monthly (-12) '
219            DO ji = 1, jpfile
220                WRITE(numout,*) trim(clname(ji)), ' frequency:', freqh(ji)
221            END DO
222         ENDIF
223
224         !--- Initialization to zero
225         qsr_oce (:,:) = 0.e0
226         qsb_oce (:,:) = 0.e0
227         qla_oce (:,:) = 0.e0
228         qlw_oce (:,:) = 0.e0
229         qnsr_oce(:,:) = 0.e0
230         qsr_ice (:,:) = 0.e0
231         qnsr_ice(:,:) = 0.e0
232         qla_ice (:,:) = 0.e0
233         qlw_ice (:,:) = 0.e0
234         qsb_ice (:,:) = 0.e0
235
236         dqns_ice(:,:) = 0.e0
237         dqla_ice(:,:) = 0.e0
238         tprecip (:,:) = 0.e0
239         sprecip (:,:) = 0.e0
240         evap    (:,:) = 0.e0
241
242         flxnow  (:,:,:) = 1.e-15         !RB
243         flxdta  (:,:,:,:) = 1.e-15       !RB
244
245         nrecflx (:)    = 0      ! switch for reading flux data for each file
246         nrecflx2(:)    = 0      ! switch for reading flux data for each file
247
248         !--- Open the files of the list
249         DO ji = 1, jpfile
250            CALL iom_open( clname(ji), numflxall(ji) )
251         END DO
252
253      ENDIF
254
255      !                                         ! ------------------------ !
256      !                                         !  Read files if required  !
257      !                                         ! ------------------------ !
258
259      !--- read data if necessary checks each file in turn
260
261      DO ji = 1, jpfile
262
263         IF (  freqh(ji) > 0 .AND.  freqh(ji) <  24 ) THEN          !--- Case of 6-hourly flux data 
264            !--- calculate current snapshot from hour of day
265            isnap = ihour / INT( freqh(ji) ) + 1
266            !--- reading is necessary when nrecflx(ji) differs from isnap
267            IF( nrecflx(ji) /= isnap ) THEN   
268               nrecflx(ji) = isnap
269               irecflx = (nday_year-1)*24 / freqh(ji) + isnap
270               irectot = 365*24 / freqh(ji)        !  this variable is not used by iom_get
271               IF(lwp) WRITE (numout,*)' CORE flux 6-h record :',irecflx, ' var:',trim(clvarname(ji))
272               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
273            ENDIF
274
275         ELSE IF ( freqh(ji) ==  24 ) THEN                             !--- Case of daily flux data 
276
277            !--- reading is necessary when nrecflx(ji) differs from nday
278            IF( nrecflx(ji) /= nday ) THEN   
279               nrecflx(ji) = nday          !! remember present read day
280               irecflx     = nday_year     !!
281               irectot = 365        !  this variable is not used by iom_get
282               IF(lwp) WRITE (numout,*)'DAILY CORE flux record :',irecflx, ' var:',trim(clvarname(ji))
283               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
284            ENDIF
285
286         ELSE IF ( freqh(ji) == -12 ) THEN                            !--- Case monthly data from CORE
287            ! we read two months all the time although we
288            ! could only read one and swap the arrays
289            IF( kt == nit000 .OR. imois /= nrecflx(ji) ) THEN
290               ! calendar computation
291               ! nrecflx number of the first file record used in the simulation
292               ! nrecflx2 number of the last  file record
293               nrecflx(ji) = imois
294               nrecflx2(ji) = nrecflx(ji)+1
295               nrecflx(ji) = MOD( nrecflx(ji), iman )
296               IF( nrecflx(ji) == 0 )   nrecflx(ji) = iman
297               nrecflx2(ji) = MOD( nrecflx2(ji), iman )
298               IF ( nrecflx2(ji) == 0 ) nrecflx2(ji) = iman
299               irectot = 12        !  this variable is not used by iom_get
300               IF(lwp) WRITE(numout,*) 'MONTHLY CORE flux record 1:', nrecflx(ji),   &
301                  &                    ' rec 2:', nrecflx2(ji), ' var:', trim(clvarname(ji))
302               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), nrecflx (ji) )
303               CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,2), nrecflx2(ji) )
304            ENDIF
305         ENDIF
306      END DO    !   end of loop over forcing files to verify if reading is necessary
307
308      !--- Printout if required     
309      !  IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
310      !            printout at the initial time step only (unless you want to debug!!)
311      IF( kt == nit000 ) THEN
312         IF(lwp) THEN
313            WRITE(numout,*)
314            WRITE(numout,*) ' read daily and monthly CORE fluxes: ok'
315            WRITE(numout,*)
316            ifpr = INT(jpi/8)      ;      jfpr = INT(jpj/10)
317            DO ji = 1, jpfile
318               WRITE(numout,*) trim(clvarname(ji)),' day: ',ndastp
319               CALL prihre(flxdta(1,1,ji,1),jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout)
320               WRITE(numout,*)
321            END DO
322         ENDIF
323      ENDIF
324
325      !                                         ! ------------------------ !
326      !                                         !   Interpolates in time   !
327      !                                         ! ------------------------ !
328      ! N.B. For now, only monthly values are interpolated, and they
329      !      are interpolated to the current day, not to the time step.
330      !
331      DO ji = 1, jpfile       
332         IF( freqh(ji) == -12 ) THEN
333            zxy = REAL(nday) / REAL( nobis(imois) ) + 0.5 - i15         !!! <== Caution nobis hard coded !!!
334            flxnow(:,:,ji) = ( ( 1. - zxy) * flxdta(:,:,ji,1) + zxy * flxdta(:,:,ji,2) )
335         ELSE
336            flxnow(:,:,ji) = flxdta(:,:,ji,1)
337         ENDIF
338      END DO
339      ! JMM : add vatm needed in tracer routines  ==> wind module ???
340      vatm(:,:) = SQRT( flxnow(:,:,2) * flxnow(:,:,2) + flxnow(:,:,3) * flxnow(:,:,3) )
341
342
343      !! =============================================== !!
344      !!   PART B:   CORE BULK CALCULATION               !!
345      !! =============================================== !!
346
347      !  for other forcing cases this is done in modules bulk.F90 and flxblk
348
349      !                                         ! ------------------------ !
350      !                                         !   Bulk initialisation    !
351      !                                         ! ------------------------ !
352
353      ! code part from bulk.F90 :
354
355      IF( kt == nit000) THEN
356         !--- computation of rdtbs2     ===>  !gm  is it really usefull ????
357         IF( nacc == 1 ) THEN
358            rdtbs2 = nfbulk * rdtmin * 0.5
359         ELSE
360           rdtbs2 = nfbulk * rdt * 0.5
361         ENDIF
362         IF( .NOT.ln_rstart ) THEN
363            zcof =  REAL( nfbulk - 1, wp )
364            gsst(:,:) =  zcof * ( tn(:,:,1)  + rt0 )
365            gsss(:,:) =  zcof * tn_ice(:,:)
366            gu  (:,:) =  zcof * un(:,:,1)
367            gv  (:,:) =  zcof * vn(:,:,1)
368         ENDIF
369      ENDIF
370
371      ! ----------------------------------------------------------------------------- !
372      !      0   Mean SST, ice temperature and ocean currents over nfbulk pdt         !
373      ! ----------------------------------------------------------------------------- !
374      !
375      gsst(:,:) = gsst(:,:) + (tn(:,:,1) + rt0 ) 
376      gsss(:,:) = gsss(:,:) + tn_ice(:,:)
377      gu  (:,:) = gu  (:,:) + un    (:,:,1)
378      gv  (:,:) = gv  (:,:) + vn    (:,:,1)
379
380      IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
381         !
382!        zsst(:,:) = gsst(:,:) / REAL( nfbulk, wp ) * tmask(:,:,1) ! mean sst in K
383!        zsss(:,:) = gsss(:,:) / REAL( nfbulk ) * tmask(:,:,1) ! mean tn_ice in K
384
385!        where ( zsst(:,:) .eq. 0.e0 ) zsst(:,:) = rt0     !lb vilain !!???
386!        where ( zsss(:,:) .eq. 0.e0 ) zsss(:,:) = rt0     !lb    //
387
388!!gm  better coded:
389         ! Caution set to rt0 over land, not 0 Kelvin (otherwise floating point exception in bulk computation)
390         zcof =  1. / REAL( nfbulk, wp )
391         zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1))    ! mean sst in K
392         zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1))   ! mean tn_ice in K
393!!gm
394
395         zut(:,:) = 0.e0      !!gm not necessary but at least first and last column
396         zvt(:,:) = 0.e0
397         !            lb
398         ! Interpolation of surface current at T-point, zut and zvt :
399!        DO ji=2, jpi
400!           zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk )
401!        END DO
402!        !
403!        DO jj=2, jpj
404!           zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk )
405!        END DO
406!!gm better coded
407         zcof =  0.5 / REAL( nfbulk, wp )
408         DO jj = 2, jpjm1
409            DO ji = 1, jpi
410               zut(ji,jj) = ( gu(ji-1,jj  ) + gu(ji,jj) ) * zcof
411               zvt(ji,jj) = ( gv(ji  ,jj-1) + gv(ji,jj) ) * zcof
412            END DO
413         END DO
414!!gm
415
416         CALL lbc_lnk( zut(:,:), 'T', -1. )         
417         CALL lbc_lnk( zvt(:,:), 'T', -1. )         
418
419         ! ----------------------------------------------------------------------------- !
420         !      I   Radiative FLUXES                                                     !
421         ! ----------------------------------------------------------------------------- !
422
423         !--- Ocean & Ice Albedo
424         alb_oce_os(:,:) = 0.   ;   alb_oce_cs(:,:) = 0.   !gm is it necessary ???
425         alb_ice_os(:,:) = 0.   ;   alb_ice_cs(:,:) = 0.
426         CALL flx_blk_albedo( alb_ice_os, alb_oce_os, alb_ice_cs, alb_oce_cs )
427
428         !--- Radiative fluxes over ocean
429         qsr_oce(:,:) = ( 1. - 0.066 ) * flxnow(:,:,5)                                 ! Solar Radiation
430         qlw_oce(:,:) = flxnow(:,:,6) - Stef*zsst(:,:)*zsst(:,:)*zsst(:,:)*zsst(:,:)   ! Long Waves (Infra-red)
431
432         !--- Radiative fluxes over ice
433         qsr_ice(:,:) = ( 1. - alb_ice_cs(:,:) ) * flxnow(:,:,5)                                  ! Solar Radiation
434         qlw_ice(:,:) = 0.95 * ( flxnow(:,:,6) - Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)*zsss(:,:) )   ! Long Waves (Infra-red)
435         !-------------------------------------------------------------------------------!
436
437         ! ----------------------------------------------------------------------------- !
438         !     II    Turbulent FLUXES                                                    !
439         ! ----------------------------------------------------------------------------- !
440         !
441         ! scalar wind ( = | U10m - SSU | )
442         ! It is important to take into account the sea surface courant
443         !                lb
444         ! Now, wind components are provided on T-points within the netcdf input file.
445         ! Thus, the wind module is computded at T-points taking into account the sea
446
447         !--- Module of relative wind
448!        dUnormt(:,:) = sqrt( (flxnow(:,:,2) - zut(:,:))*(flxnow(:,:,2) - zut(:,:))  &
449!             &             + (flxnow(:,:,3) - zvt(:,:))*(flxnow(:,:,3) - zvt(:,:)) )
450!!gm  more efficient coding:
451         DO jj = 1, jpj
452            DO ji = 1, jpi
453               zzu = flxnow(ji,jj,2) - zut(ji,jj)
454               zzv = flxnow(ji,jj,3) - zvt(ji,jj)
455               dUnormt(ji,jj) = SQRT( zzu*zzu + zzv*zzv )
456            END DO
457         END DO
458!!gm
459         !--- specific humidity at temp SST
460         qsatw(:,:) = 0.98*640380*exp(-5107.4/zsst(:,:))/rhoa
461         !
462         !--- specific humidity at temp tn_ice
463         qsat(:,:)  = 11637800*exp(-5897.8/zsss(:,:))/rhoa
464         !
465
466         ! CORE iterartive algo for computation of Cd, Ch, Ce at T-point :
467         ! ===============================================================
468         IF( ln_2m ) THEN
469            IF( kt == nit000 .AND. lwp )   THEN
470               WRITE(numout,*)
471               WRITE(numout,*) '           Calling TURB_CORE_2Z for bulk transfert coefficients'
472            ENDIF
473            !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) :
474            CALL TURB_CORE_2Z(2.,10.,zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), &
475               &              dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:), zt_zu(:,:), zq_zu(:,:))
476         ELSE
477            IF( kt == nit000 .AND. lwp )    THEN
478               WRITE(numout,*)
479               WRITE(numout,*) '           Calling TURB_CORE_1Z for bulk transfert coefficients'
480            ENDIF
481            !! If air temp. and spec. hum. are given at same height than wind (10m) :
482            CALL TURB_CORE_1Z(10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4),  &
483               &              dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) )
484         ENDIF
485
486         !   II.1) Momentum over ocean and ice
487         !   ---------------------------------
488         ! Tau_x at T-point
489         tauxt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,2) - zut(:,:))    &
490              &                          + freeze(:,:)*Cice*flxnow(:,:,2) ) !lb correct pour glace
491         ! Tau_y at T-point
492         tauyt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,3) - zvt(:,:))    &
493              &                          + freeze(:,:)*Cice*flxnow(:,:,3) )
494         !
495         CALL lbc_lnk( tauxt(:,:), 'T', -1. )       !!gm seems not decessary at this point....
496         CALL lbc_lnk( tauyt(:,:), 'T', -1. )
497
498         !
499         !Tau_x at U-point
500         DO ji = 1, jpim1
501            taux(ji,:) = 0.5*(tauxt(ji,:) + tauxt(ji+1,:))
502         END DO
503         !
504         ! Tau_y at V-point
505         DO jj = 1, jpjm1
506            tauy(:,jj) = 0.5*(tauyt(:,jj) + tauyt(:,jj+1))
507         END DO
508
509         !
510         !   II.2) Turbulent fluxes over ocean
511         !   ---------------------------------
512         !
513         IF( ln_2m ) THEN
514            !!
515            !! Values of temp. and hum. adjusted to 10m must be used instead of 2m values
516            !! Sensible Heat :    ! right sign for ocean
517            qsb_oce(:,:) = rhoa*cp*Ch(:,:)*(zt_zu(:,:) - zsst(:,:))*dUnormt(:,:)
518            !!
519            !! Latent Heat :      ! wrong sign for ocean
520            evap(:,:)    = rhoa*Ce(:,:)*(qsatw(:,:) - zq_zu(:,:))*dUnormt(:,:)
521            !!
522         ELSE
523            !!
524            !! Sensible Heat :    ! right sign for ocean
525            qsb_oce(:,:) = rhoa*cp*Ch(:,:)*(flxnow(:,:,7) - zsst(:,:))*dUnormt(:,:)
526            !!
527            !! Latent Heat :      ! wrong sign for ocean
528            evap(:,:)    = rhoa*Ce(:,:)*(qsatw(:,:) - flxnow(:,:,4))*dUnormt(:,:)
529            !!
530         END IF
531         !!
532         qla_oce(:,:) =  -Lv * evap(:,:) ! right sign for ocean
533         
534         !
535         !   II.3) Turbulent fluxes over ice
536         !   -------------------------------
537         !
538         ! Sensible Heat :
539         qsb_ice(:,:) = rhoa*cp*Cice*( flxnow(:,:,7) - zsss(:,:) )*dUnormt(:,:) !lb     use
540         !                                                                    !lb   dUnormt???
541         ! Latent Heat :                                                      !lb or rather Unormt?
542         qla_ice(:,:) = Ls*rhoa*Cice*( flxnow(:,:,4) - qsat(:,:) )*dUnormt(:,:)
543         !
544         !-------------------------------------------------------------------------------------
545
546
547         ! ----------------------------------------------------------------------------- !
548         !     III    Total FLUXES                                                       !
549         ! ----------------------------------------------------------------------------- !
550         !
551         !   III.1) Downward Non Solar flux over ocean
552         !   -----------------------------------------
553         qnsr_oce(:,:) = qlw_oce(:,:) + qsb_oce(:,:) + qla_oce(:,:)
554         !
555         !   III.1) Downward Non Solar flux over ice
556         !   ---------------------------------------
557         qnsr_ice(:,:) = qlw_ice(:,:) + qsb_ice(:,:) + qla_ice(:,:)
558         !
559         !-------------------------------------------------------------------------------!
560
561         ! 6. TOTAL NON SOLAR SENSITIVITY
562
563         dqlw_ice(:,:)= 4.0*0.95*Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)
564
565         ! d qla_ice/ d zsss
566         dqla_ice(:,:) = -Ls*Cice*11637800 &
567              * (-5897.8)/(zsss(:,:)*zsss(:,:)) &
568              * exp(-5897.8/zsss(:,:)) * dUnormt(:,:)
569
570         ! d qsb_ice/ d zsss
571         dqsb_ice(:,:) = rhoa * cp * Cice * dUnormt(:,:)
572
573         dqns_ice(:,:) = - ( dqlw_ice(:,:) + dqsb_ice(:,:) + dqla_ice(:,:) )
574
575         !--------------------------------------------------------------------
576         ! FRACTION of net shortwave radiation which is not absorbed in the
577         ! thin surface layer and penetrates inside the ice cover
578         ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
579
580         catm1(:,:) = 1.0 - 0.3
581         fr1_i0(:,:) = 0.18 * ( 1.0 - catm1(:,:) ) + 0.35 * catm1(:,:)
582         fr2_i0(:,:) = 0.82 * ( 1.0 - catm1(:,:) ) + 0.65 * catm1(:,:)
583
584         !  Total  PRECIPITATION (kg/m2/s)
585         tprecip(:,:) = alpha_precip*flxnow(:,:,1)
586         ! rename precipitation for freshwater budget calculations:
587!         watm(:,:) =  flxnow(:,:,1)*1000
588         watm(:,:) =  flxnow(:,:,1)*rday
589         !
590
591         !  SNOW PRECIPITATION  (kg/m2/s)
592         sprecip(:,:) =  alpha_precip*flxnow(:,:,8)
593
594         !---------------------------------------------------------------------
595
596         CALL lbc_lnk( taux (:,:)    , 'U', -1. )
597         CALL lbc_lnk( tauy (:,:)    , 'V', -1. )
598         CALL lbc_lnk( qsr_oce (:,:) , 'T', 1. )
599         CALL lbc_lnk( qnsr_oce(:,:) , 'T', 1. )
600         CALL lbc_lnk( qsr_ice (:,:) , 'T', 1. )
601         CALL lbc_lnk( qnsr_ice(:,:) , 'T', 1. )
602         CALL lbc_lnk( qla_ice (:,:) , 'T', 1. )
603         CALL lbc_lnk( dqns_ice(:,:) , 'T', 1. )
604         CALL lbc_lnk( dqla_ice(:,:) , 'T', 1. )
605         CALL lbc_lnk( fr1_i0  (:,:) , 'T', 1. )
606         CALL lbc_lnk( fr2_i0  (:,:) , 'T', 1. )
607         CALL lbc_lnk( tprecip (:,:) , 'T', 1. )
608         CALL lbc_lnk( sprecip (:,:) , 'T', 1. )
609         CALL lbc_lnk( evap    (:,:) , 'T', 1. )
610         !!
611         !!
612         !!   NEVER mask the windstress!!
613         qsr_oce (:,:) = qsr_oce (:,:)*tmask(:,:,1)
614         qnsr_oce(:,:) = qnsr_oce(:,:)*tmask(:,:,1)
615         qsr_ice (:,:) = qsr_ice (:,:)*tmask(:,:,1)
616         qnsr_ice(:,:) = qnsr_ice(:,:)*tmask(:,:,1)
617         qla_ice (:,:) = qla_ice (:,:)*tmask(:,:,1)
618         dqns_ice(:,:) = dqns_ice(:,:)*tmask(:,:,1)
619         dqla_ice(:,:) = dqla_ice(:,:)*tmask(:,:,1)
620         fr1_i0  (:,:) = fr1_i0  (:,:)*tmask(:,:,1)
621         fr2_i0  (:,:) = fr2_i0  (:,:)*tmask(:,:,1)
622         tprecip (:,:) = tprecip (:,:)*tmask(:,:,1)
623         sprecip (:,:) = sprecip (:,:)*tmask(:,:,1)
624         evap    (:,:) = evap    (:,:)*tmask(:,:,1)
625         gsst(:,:) = 0.e0
626         gsss(:,:) = 0.e0
627         gu  (:,:) = 0.e0
628         gv  (:,:) = 0.e0
629
630         !   Degugging print (A.M. Treguier)       
631         !      write (55) taux, tauy, qnsr_oce, qsb_ice, qnsr_ice, qla_ice, dqns_ice &
632         !  &     , dqla_ice, fr1_i0, fr2_i0, qlw_oce, qla_oce, qsb_oce, evap
633         !      write(numout,*) 'written 14 arrays on unit fort.55'
634
635
636      ENDIF
637
638      ! ------------------- !
639      ! Last call kt=nitend !
640      ! ------------------- !
641
642      ! Closing of the numflx file (required in mpp)
643
644      IF( kt == nitend ) THEN
645         DO ji=1, jpfile
646            CALL iom_close(numflxall(ji))
647         END DO
648      ENDIF
649      IF ( lrst_oce ) CALL core_rst( kt, 'WRITE' )
650
651   END SUBROUTINE flx
652
653   SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   &
654      &                        dU, Cd, Ch, Ce   )
655      !!----------------------------------------------------------------------
656      !!                      ***  ROUTINE  turb_core  ***
657      !!
658      !! ** Purpose :   Computes turbulent transfert coefficients of surface
659      !!                fluxes according to Large & Yeager (2004)
660      !!
661      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
662      !!      Momentum, Latent and sensible heat exchange coefficients
663      !!      Caution: this procedure should only be used in cases when air
664      !!      temperature (T_air), air specific humidity (q_air) and wind (dU)
665      !!      are provided at the same height 'zzu'!
666      !!
667      !! References :
668      !!      Large & Yeager, 2004 : ???
669      !! History :
670      !!        !  XX-XX  (???     )  Original code
671      !!   9.0  !  05-08  (L. Brodeau) Rewriting and optimization
672      !!----------------------------------------------------------------------
673      !! * Arguments
674
675      REAL(wp), INTENT(in) :: zu                 ! altitude of wind measurement       [m]
676      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::  &
677         sst,       &       ! sea surface temperature         [Kelvin]
678         T_a,       &       ! potential air temperature       [Kelvin]
679         q_sat,     &       ! sea surface specific humidity   [kg/kg]
680         q_a,       &       ! specific air humidity           [kg/kg]
681         dU                 ! wind module |U(zu)-U(0)|        [m/s]
682      REAL(wp), intent(out), DIMENSION(jpi,jpj)  :: &
683         Cd,    &                ! transfert coefficient for momentum       (tau)
684         Ch,    &                ! transfert coefficient for temperature (Q_sens)
685         Ce                      ! transfert coefficient for evaporation  (Q_lat)
686
687      !! * Local declarations
688      REAL(wp), DIMENSION(jpi,jpj)  ::   &
689         dU10,        &       ! dU                                   [m/s]
690         dT,          &       ! air/sea temperature differeence      [K]
691         dq,          &       ! air/sea humidity difference          [K]
692         Cd_n10,      &       ! 10m neutral drag coefficient
693         Ce_n10,      &       ! 10m neutral latent coefficient
694         Ch_n10,      &       ! 10m neutral sensible coefficient
695         sqrt_Cd_n10, &       ! root square of Cd_n10
696         sqrt_Cd,     &       ! root square of Cd
697         T_vpot,      &       ! virtual potential temperature        [K]
698         T_star,      &       ! turbulent scale of tem. fluct.
699         q_star,      &       ! turbulent humidity of temp. fluct.
700         U_star,      &       ! turb. scale of velocity fluct.
701         L,           &       ! Monin-Obukov length                  [m]
702         zeta,        &       ! stability parameter at height zu
703         U_n10,       &       ! neutral wind velocity at 10m         [m]   
704         xlogt, xct, zpsi_h, zpsi_m
705      !!
706      INTEGER :: j_itt
707      INTEGER, PARAMETER :: nb_itt = 3
708      INTEGER, DIMENSION(jpi,jpj)  ::   &
709         stab         ! 1st guess stability test integer
710
711      REAL(wp), PARAMETER ::                        &
712         grav   = 9.8,          &  ! gravity                       
713         kappa  = 0.4              ! von Karman s constant
714
715      !! * Start
716      !! Air/sea differences
717      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s
718      dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu
719      dq = q_a - q_sat
720      !!   
721      !! Virtual potential temperature
722      T_vpot = T_a*(1. + 0.608*q_a)
723      !!
724      !! Neutral Drag Coefficient
725      stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0
726      Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
727      sqrt_Cd_n10 = sqrt(Cd_n10)
728      Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b)
729      Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !   L & Y eq. (6c), (6d)
730      !!
731      !! Initializing transfert coefficients with their first guess neutral equivalents :
732      Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd)
733
734      !! * Now starting iteration loop
735      DO j_itt=1, nb_itt
736         !! Turbulent scales :
737         U_star  = sqrt_Cd*dU10                !   L & Y eq. (7a)
738         T_star  = Ch/sqrt_Cd*dT               !   L & Y eq. (7b)
739         q_star  = Ce/sqrt_Cd*dq               !   L & Y eq. (7c)
740
741         !! Estimate the Monin-Obukov length :
742         L  = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) )
743
744         !! Stability parameters :
745         zeta = zu/L ;  zeta   = sign( min(abs(zeta),10.0), zeta )
746         zpsi_h = psi_h(zeta)
747         zpsi_m = psi_m(zeta)
748
749         !! Shifting the wind speed to 10m and neutral stability :
750         U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a)
751
752         !! Updating the neutral 10m transfer coefficients :
753         Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
754         sqrt_Cd_n10 = sqrt(Cd_n10)
755         Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
756         stab    = 0.5 + sign(0.5,zeta)
757         Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d)
758
759         !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
760         !!
761         xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m)
762         Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd)
763         !!
764         xlogt = log(zu/10.) - zpsi_h
765         !!
766         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
767         Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
768         !!
769         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
770         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
771         !!
772      END DO
773      !!
774    END SUBROUTINE TURB_CORE_1Z
775
776
777    SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu)
778      !!----------------------------------------------------------------------
779      !!                      ***  ROUTINE  turb_core  ***
780      !!
781      !! ** Purpose :   Computes turbulent transfert coefficients of surface
782      !!                fluxes according to Large & Yeager (2004).
783      !!
784      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
785      !!      Momentum, Latent and sensible heat exchange coefficients
786      !!      Caution: this procedure should only be used in cases when air
787      !!      temperature (T_air) and air specific humidity (q_air) are at 2m
788      !!      whereas wind (dU) is at 10m.
789      !!
790      !! References :
791      !!      Large & Yeager, 2004 : ???
792      !! History :
793      !!   9.0  !  06-12  (L. Brodeau) Original code for 2Z
794      !!----------------------------------------------------------------------
795      !! * Arguments
796      REAL(wp), INTENT(in)   :: &
797         zt,      &     ! height for T_zt and q_zt                   [m]
798         zu             ! height for dU                              [m]
799      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::  &
800         sst,      &     ! sea surface temperature              [Kelvin]
801         T_zt,     &     ! potential air temperature            [Kelvin]
802         q_sat,    &     ! sea surface specific humidity         [kg/kg]
803         q_zt,     &     ! specific air humidity                 [kg/kg]
804         dU              ! relative wind module |U(zu)-U(0)|       [m/s]
805      REAL(wp), INTENT(out), DIMENSION(jpi,jpj)  ::  &
806         Cd,       &     ! transfer coefficient for momentum         (tau)
807         Ch,       &     ! transfer coefficient for sensible heat (Q_sens)
808         Ce,       &     ! transfert coefficient for evaporation   (Q_lat)
809         T_zu,     &     ! air temp. shifted at zu                     [K]
810         q_zu            ! spec. hum.  shifted at zu               [kg/kg]
811
812      !! * Local declarations
813      REAL(wp), DIMENSION(jpi,jpj) ::  &
814         dU10,        &     ! dU                                [m/s]
815         dT,          &     ! air/sea temperature differeence   [K]
816         dq,          &     ! air/sea humidity difference       [K]
817         Cd_n10,      &     ! 10m neutral drag coefficient
818         Ce_n10,      &     ! 10m neutral latent coefficient
819         Ch_n10,      &     ! 10m neutral sensible coefficient
820         sqrt_Cd_n10, &     ! root square of Cd_n10
821         sqrt_Cd,     &     ! root square of Cd
822         T_vpot_u,    &     ! virtual potential temperature        [K]
823         T_star,      &     ! turbulent scale of tem. fluct.
824         q_star,      &     ! turbulent humidity of temp. fluct.
825         U_star,      &     ! turb. scale of velocity fluct.
826         L,           &     ! Monin-Obukov length                  [m]
827         zeta_u,      &     ! stability parameter at height zu
828         zeta_t,      &     ! stability parameter at height zt
829         U_n10,       &     ! neutral wind velocity at 10m        [m]
830         xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m
831
832      INTEGER :: j_itt
833      INTEGER, PARAMETER :: nb_itt = 3   ! number of itterations
834      INTEGER, DIMENSION(jpi,jpj) :: &
835           &     stab                ! 1st stability test integer
836      REAL(wp), PARAMETER ::                        &
837         grav   = 9.8,      &  ! gravity                       
838         kappa  = 0.4          ! von Karman's constant
839
840      !!  * Start
841
842      !! Initial air/sea differences
843      dU10 = max(0.5, dU)      !  we don't want to fall under 0.5 m/s
844      dT = T_zt - sst
845      dq = q_zt - q_sat
846
847      !! Neutral Drag Coefficient :
848      stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE
849      Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 )
850      sqrt_Cd_n10 = sqrt(Cd_n10)
851      Ce_n10  = 1E-3*( 34.6 * sqrt_Cd_n10 )
852      Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab))
853
854      !! Initializing transf. coeff. with their first guess neutral equivalents :
855      Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd)
856
857      !! Initializing z_u values with z_t values :
858      T_zu = T_zt ;  q_zu = q_zt
859
860      !!  * Now starting iteration loop
861      DO j_itt=1, nb_itt
862         dT = T_zu - sst ;  dq = q_zu - q_sat ! Updating air/sea differences
863         T_vpot_u = T_zu*(1. + 0.608*q_zu)    ! Updating virtual potential temperature at zu
864         U_star = sqrt_Cd*dU10                ! Updating turbulent scales :   (L & Y eq. (7))
865         T_star  = Ch/sqrt_Cd*dT              !
866         q_star  = Ce/sqrt_Cd*dq              !
867         !!
868         L = (U_star*U_star) &                ! Estimate the Monin-Obukov length at height zu
869              & / (kappa*grav/T_vpot_u*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star))
870         !! Stability parameters :
871         zeta_u  = zu/L  ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
872         zeta_t  = zt/L  ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t )
873         zpsi_hu = psi_h(zeta_u)
874         zpsi_ht = psi_h(zeta_t)
875         zpsi_m  = psi_m(zeta_u)
876         !!
877         !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a))
878!        U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)))
879         U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m))
880         !!
881         !! Shifting temperature and humidity at zu :          (L & Y eq. (9b-9c))
882!        T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
883         T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
884!        q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
885         q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
886         !!
887         !! q_zu cannot have a negative value : forcing 0
888         stab = 0.5 + sign(0.5,q_zu) ;  q_zu = stab*q_zu
889         !!
890         !! Updating the neutral 10m transfer coefficients :
891         Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
892         sqrt_Cd_n10 = sqrt(Cd_n10)
893         Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
894         stab    = 0.5 + sign(0.5,zeta_u)
895         Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d)
896         !!
897         !!
898         !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
899!        xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))
900         xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)
901         Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
902         !!
903!        xlogt = log(zu/10.) - psi_h(zeta_u)
904         xlogt = log(zu/10.) - zpsi_hu
905         !!
906         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
907         Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
908         !!
909         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
910         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
911         !!
912         !!
913      END DO
914      !!
915    END SUBROUTINE TURB_CORE_2Z
916
917    FUNCTION psi_m(zta)   !! Psis, L & Y eq. (8c), (8d), (8e)
918      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
919
920      REAL(wp), PARAMETER :: pi = 3.14159
921      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m
922      REAL(wp), DIMENSION(jpi,jpj)             :: X2, X, stabit
923      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.0) ;  X  = sqrt(X2)
924      stabit    = 0.5 + sign(0.5,zta)
925      psi_m = -5.*zta*stabit  &                                                  ! Stable
926           & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2)  ! Unstable
927    END FUNCTION psi_m
928
929    FUNCTION psi_h(zta)    !! Psis, L & Y eq. (8c), (8d), (8e)
930      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
931
932      REAL(wp), DIMENSION(jpi,jpj)             :: psi_h
933      REAL(wp), DIMENSION(jpi,jpj)             :: X2, X, stabit
934      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.) ;  X  = sqrt(X2)
935      stabit    = 0.5 + sign(0.5,zta)
936      psi_h = -5.*zta*stabit  &                                       ! Stable
937           & + (1. - stabit)*(2.*log( (1. + X2)/2. ))                 ! Unstable
938    END FUNCTION psi_h
939
940
941  SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp )
942    !!----------------------------------------------------------------------
943    !!               ***  ROUTINE flx_blk_albedo  ***
944    !!
945    !! ** Purpose :   Computation of the albedo of the snow/ice system
946    !!      as well as the ocean one
947    !!
948    !! ** Method  : - Computation of the albedo of snow or ice (choose the
949    !!      right one by a large number of tests
950    !!              - Computation of the albedo of the ocean
951    !!
952    !! References :
953    !!      Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
954    !!
955    !! History :
956    !!  8.0   !  01-04  (LIM 1.0)
957    !!  8.5   !  03-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
958    !!----------------------------------------------------------------------
959    !! * Modules used
960    USE ice                   ! ???
961
962    !! * Arguments
963    REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::  &
964         palb         , & !  albedo of ice under overcast sky
965         palcn        , & !  albedo of ocean under overcast sky
966         palbp        , & !  albedo of ice under clear sky
967         palcnp                  !  albedo of ocean under clear sky
968
969    !! * Local variables
970    INTEGER ::    &
971         ji, jj                   ! dummy loop indices
972    REAL(wp) ::   &
973         c1     = 0.05  , & ! constants values
974         c2     = 0.1   , &
975         albice = 0.50  , & !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
976         cgren  = 0.06  , & !  correction of the snow or ice albedo to take into account
977         !  effects of cloudiness (Grenfell & Perovich, 1984)
978         alphd  = 0.80  , & !  coefficients for linear interpolation used to compute
979         alphdi = 0.72  , & !  albedo between two extremes values (Pyane, 1972)
980         alphc  = 0.65  , &
981         zmue   = 0.4   , &     !  cosine of local solar altitude
982         zzero   = 0.0     ,  &
983         zone    = 1.0
984
985    REAL(wp) ::   &
986         zmue14         , & !  zmue**1.4
987         zalbpsnm       , & !  albedo of ice under clear sky when snow is melting
988         zalbpsnf       , & !  albedo of ice under clear sky when snow is freezing
989         zalbpsn        , & !  albedo of snow/ice system when ice is coverd by snow
990         zalbpic        , & !  albedo of snow/ice system when ice is free of snow
991         zithsn         , & !  = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
992         zitmlsn        , & !  = 1 freezinz snow (sist >=rt0_snow) ; = 0 melting snow (sist<rt0_snow)
993         zihsc1         , & !  = 1 hsn <= c1 ; = 0 hsn > c1
994         zihsc2                   !  = 1 hsn >= c2 ; = 0 hsn < c2
995    REAL(wp), DIMENSION(jpi,jpj) ::  &
996         zalbfz         , & !  ( = alphdi for freezing ice ; = albice for melting ice )
997         zficeth                  !  function of ice thickness
998    LOGICAL , DIMENSION(jpi,jpj) ::  &
999         llmask
1000    !!    to be included for without seaice
1001    !!      REAL(wp), DIMENSION(jpi,jpj) ::   &  !:
1002    !!      sist   ,   &  !: Sea-Ice Surface Temperature (Kelvin )
1003    !!      hsnif  ,   &  !: Snow thickness
1004    !!      hicif       !: Ice thickness
1005   
1006    !!---------------------------------------------------------------------
1007
1008    !-------------------------
1009    !  Computation of  zficeth
1010    !--------------------------
1011
1012    llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice )
1013    WHERE ( llmask )   !  ice free of snow and melts
1014       zalbfz = albice
1015    ELSEWHERE
1016       zalbfz = alphdi
1017    END WHERE
1018
1019    DO jj = 1, jpj
1020       DO ji = 1, jpi
1021          IF( hicif(ji,jj) > 1.5 ) THEN
1022             zficeth(ji,jj) = zalbfz(ji,jj)
1023          ELSEIF( hicif(ji,jj) > 1.0  .AND. hicif(ji,jj) <= 1.5 ) THEN
1024             zficeth(ji,jj) = 0.472 + 2.0 * ( zalbfz(ji,jj) - 0.472 ) * ( hicif(ji,jj) - 1.0 )
1025          ELSEIF( hicif(ji,jj) > 0.05 .AND. hicif(ji,jj) <= 1.0 ) THEN
1026             zficeth(ji,jj) = 0.2467 + 0.7049 * hicif(ji,jj)                                &
1027                  &                    - 0.8608 * hicif(ji,jj) * hicif(ji,jj)                 &
1028                  &                    + 0.3812 * hicif(ji,jj) * hicif(ji,jj) * hicif (ji,jj)
1029          ELSE
1030             zficeth(ji,jj) = 0.1 + 3.6 * hicif(ji,jj)
1031          ENDIF
1032       END DO
1033    END DO
1034
1035    !-----------------------------------------------
1036    !    Computation of the snow/ice albedo system
1037    !-------------------------- ---------------------
1038
1039    !    Albedo of snow-ice for clear sky.
1040    !-----------------------------------------------
1041     DO jj = 1, jpj
1042        DO ji = 1, jpi
1043           !  Case of ice covered by snow.
1044           !  melting snow
1045           zihsc1       = 1.0 - MAX ( zzero , SIGN ( zone , - ( hsnif(ji,jj) - c1 ) ) )
1046           zalbpsnm     = ( 1.0 - zihsc1 ) &
1047                * ( zficeth(ji,jj) + hsnif(ji,jj) * ( alphd - zficeth(ji,jj) ) / c1 ) &
1048                &                 + zihsc1   * alphd
1049          !  freezing snow
1050           zihsc2       = MAX ( zzero , SIGN ( zone , hsnif(ji,jj) - c2 ) )
1051           zalbpsnf     = ( 1.0 - zihsc2 ) * &
1052                ( albice + hsnif(ji,jj) * ( alphc - albice ) / c2 ) &
1053                &                 + zihsc2   * alphc
1054
1055           zitmlsn      =  MAX ( zzero , SIGN ( zone , sist(ji,jj) - rt0_snow ) )
1056           zalbpsn      =  zitmlsn * zalbpsnf + ( 1.0 - zitmlsn ) * zalbpsnm
1057
1058           !  Case of ice free of snow.
1059           zalbpic      = zficeth(ji,jj)
1060
1061           ! albedo of the system
1062           zithsn       = 1.0 - MAX ( zzero , SIGN ( zone , - hsnif(ji,jj) ) )
1063           palbp(ji,jj) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
1064        END DO
1065     END DO
1066
1067     !    Albedo of snow-ice for overcast sky.
1068     !----------------------------------------------
1069     palb(:,:)   = palbp(:,:) + cgren
1070
1071     !--------------------------------------------
1072     !    Computation of the albedo of the ocean
1073     !-------------------------- -----------------
1074
1075
1076     !  Parameterization of Briegled and Ramanathan, 1982
1077     zmue14      = zmue**1.4
1078     palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 )
1079
1080     !  Parameterization of Kondratyev, 1969 and Payne, 1972
1081     palcn(:,:)  = 0.06
1082
1083   END SUBROUTINE flx_blk_albedo
1084
1085   SUBROUTINE core_rst( kt, cdrw )
1086     !!---------------------------------------------------------------------
1087     !!                   ***  ROUTINE ts_rst  ***
1088     !!
1089     !! ** Purpose : Read or write CORE specific variables ( gsss, gu, gv )
1090     !!
1091     !! ** Method  :
1092     !!
1093     !!----------------------------------------------------------------------
1094     USE iom          ! move to top of flxmod to avoid duplication
1095     USE restart      ! move to top of flxmod to avoid duplication
1096     !
1097     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1098     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1099     !
1100     !!----------------------------------------------------------------------
1101     !
1102     IF( TRIM(cdrw) == 'READ' ) THEN
1103        IF( ln_rstart ) THEN
1104              CALL iom_get( numror, jpdom_autoglo, 'gsss', gsss )
1105              CALL iom_get( numror, jpdom_autoglo, 'gu', gu )
1106              CALL iom_get( numror, jpdom_autoglo, 'gv', gv )
1107        ENDIF ! gsss, gu, gv are initialized in the routine ( may be changed )
1108     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1109        CALL iom_rstput( kt, nitrst, numrow, 'gsss', gsss )
1110        CALL iom_rstput( kt, nitrst, numrow, 'gu', gu )
1111        CALL iom_rstput( kt, nitrst, numrow, 'gv', gv )
1112     ENDIF
1113     !
1114   END SUBROUTINE core_rst
Note: See TracBrowser for help on using the repository browser.