source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/surf_ocean_mod.F90 @ 3699

Last change on this file since 3699 was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 11.9 KB
Line 
1!
2! $Id: surf_ocean_mod.F90 2538 2016-06-03 14:12:16Z fairhead $
3!
4MODULE surf_ocean_mod
5
6  IMPLICIT NONE
7
8CONTAINS
9!
10!******************************************************************************
11!
12  SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
13       windsp, rmu0, fder, tsurf_in, &
14       itime, dtime, jour, knon, knindex, &
15       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
16       AcoefH, AcoefQ, BcoefH, BcoefQ, &
17       AcoefU, AcoefV, BcoefU, BcoefV, &
18       ps, u1, v1, gustiness, rugoro, pctsrf, &
19       snow, qsurf, agesno, &
20       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &
21       tsurf_new, dflux_s, dflux_l, lmt_bils, &
22       flux_u1, flux_v1 &
23#ifdef ISO
24        &       ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
25        &       xtsnow,xtevap,h1 &
26#endif               
27        &       )
28
29  use albedo, only: alboc, alboc_cd
30  USE dimphy, ONLY: klon, zmasq
31  USE surface_data, ONLY     : type_ocean
32  USE ocean_forced_mod, ONLY : ocean_forced_noice
33  USE ocean_slab_mod, ONLY   : ocean_slab_noice
34  USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
35  USE indice_sol_mod, ONLY : nbsrf, is_oce
36#ifdef ISO
37  USE infotrac_phy, ONLY : ntraciso,niso
38#ifdef ISOVERIF
39    USE isotopes_mod, ONLY: iso_eau,ridicule
40    USE isotopes_verif_mod
41#endif
42#endif
43!
44! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force,
45! slab or couple). The calculations of albedo and rugosity for the ocean surface are
46! done in here because they are identical for the different modes of ocean.
47
48
49    INCLUDE "YOMCST.h"
50
51    include "clesphys.h"
52    ! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0)
53
54! Input variables
55!******************************************************************************
56    INTEGER, INTENT(IN)                      :: itime, jour, knon
57    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
58    REAL, INTENT(IN)                         :: dtime
59    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
60    REAL, DIMENSION(klon), INTENT(IN)        :: swnet  ! net shortwave radiation at surface 
61    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface 
62    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
63    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
64    REAL, DIMENSION(klon), INTENT(IN)        :: rmu0 
65    REAL, DIMENSION(klon), INTENT(IN)        :: fder
66    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
67    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
68    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
69    REAL, DIMENSION(klon), INTENT(IN)        :: cdragm
70    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
71    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
72    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
73    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
74    REAL, DIMENSION(klon), INTENT(IN)        :: ps
75    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
76    REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
77    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
78#ifdef ISO
79    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: xtprecip_rain, xtprecip_snow
80    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: xtspechum
81#endif
82
83! In/Output variables
84!******************************************************************************
85    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
86    REAL, DIMENSION(klon), INTENT(INOUT)     :: qsurf
87    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
88#ifdef ISO
89    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsnow
90    REAL, DIMENSION(niso,klon), INTENT(INOUT)        :: Roce 
91#endif
92
93! Output variables
94!******************************************************************************
95    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
96!albedo SB >>>
97!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
98!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
99    REAL, DIMENSION(6), INTENT(IN)          :: SFRWL
100    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
101!albedo SB <<<     
102    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
103    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
104    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
105    REAL, DIMENSION(klon), INTENT(OUT)       :: lmt_bils
106    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
107#ifdef ISO
108    REAL, DIMENSION(ntraciso,klon), INTENT(out)        :: xtevap ! isotopes in surface evaporation flux
109    REAL, DIMENSION(klon), INTENT(out)        :: h1 ! just a diagnostic, not useful for the simulation   
110#endif
111
112! Local variables
113!******************************************************************************
114    INTEGER               :: i, k
115    REAL                  :: tmp
116    REAL, PARAMETER       :: cepdu2=(0.1)**2
117    REAL, DIMENSION(klon) :: alb_eau
118    REAL, DIMENSION(klon) :: radsol
119#ifdef ISO
120#ifdef ISOVERIF   
121    integer j
122#endif   
123#endif
124    REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
125    CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
126
127! End definition
128!******************************************************************************
129
130#ifdef ISO
131#ifdef ISOVERIF
132        do i=1,knon
133          if (iso_eau.gt.0) then         
134                 call iso_verif_egalite_choix(xtspechum(iso_eau,i), &
135     &                  spechum(i),'surf_ocean_mod 117', &
136     &                  errmax,errmaxrel)         
137                 call iso_verif_egalite_choix(xtsnow(iso_eau,i), &
138     &                  snow(i),'surf_ocean_mod 127', &
139     &                  errmax,errmaxrel)
140           endif !if (iso_eau.gt.0) then
141        enddo !do i=1,klon
142#endif     
143#endif
144
145!******************************************************************************
146! Calculate total net radiance at surface
147!
148!******************************************************************************
149    radsol(:) = 0.0
150    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
151
152!******************************************************************************
153! Cdragq computed from cdrag
154! The difference comes only from a factor (f_z0qh_oce) on z0, so that
155! it can be computed inside surf_ocean
156! More complicated appraches may require the propagation through
157! pbl_surface of an independant cdragq variable.
158!******************************************************************************
159
160    IF ( f_z0qh_oce .ne. 1.) THEN
161! Si on suit les formulations par exemple de Tessel, on
162! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
163       cdragq(:)=cdragh(:)*                                      &
164       log(z1lay(:)/z0h(:))/log(z1lay(:)/(f_z0qh_oce*z0h(:)))
165    ELSE
166       cdragq(:)=cdragh(:)
167    ENDIF
168
169!******************************************************************************
170! Switch according to type of ocean (couple, slab or forced)
171!******************************************************************************
172    SELECT CASE(type_ocean)
173    CASE('couple')
174       CALL ocean_cpl_noice( &
175            swnet, lwnet, alb1, &
176            windsp, fder, &
177            itime, dtime, knon, knindex, &
178            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,&
179            AcoefH, AcoefQ, BcoefH, BcoefQ, &
180            AcoefU, AcoefV, BcoefU, BcoefV, &
181            ps, u1, v1, gustiness, &
182            radsol, snow, agesno, &
183            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
184            tsurf_new, dflux_s, dflux_l)
185
186    CASE('slab')
187       CALL ocean_slab_noice( &
188            itime, dtime, jour, knon, knindex, &
189            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
190            AcoefH, AcoefQ, BcoefH, BcoefQ, &
191            AcoefU, AcoefV, BcoefU, BcoefV, &
192            ps, u1, v1, gustiness, tsurf_in, &
193            radsol, snow, &
194            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
195            tsurf_new, dflux_s, dflux_l, lmt_bils)
196       
197    CASE('force')
198       CALL ocean_forced_noice( &
199            itime, dtime, jour, knon, knindex, &
200            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
201            temp_air, spechum, &
202            AcoefH, AcoefQ, BcoefH, BcoefQ, &
203            AcoefU, AcoefV, BcoefU, BcoefV, &
204            ps, u1, v1, gustiness, &
205            radsol, snow, agesno, &
206            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
207            tsurf_new, dflux_s, dflux_l &
208#ifdef ISO
209            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
210            xtsnow,xtevap,h1 & 
211#endif           
212            )
213    END SELECT
214
215!******************************************************************************
216! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
217!******************************************************************************
218    IF (type_ocean.NE.'slab') THEN
219        lmt_bils(:)=0.
220        DO i=1,knon
221           lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
222           *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i)))
223        END DO
224    END IF
225
226!******************************************************************************
227! Calculate ocean surface albedo
228!******************************************************************************
229!albedo SB >>>
230IF (iflag_albedo==0) THEN
231!--old parametrizations of ocean surface albedo
232!
233    IF (cycle_diurne) THEN
234!
235       CALL alboc_cd(rmu0,alb_eau)
236!
237!--ad-hoc correction for model radiative balance tuning
238!--now outside alboc_cd routine
239       alb_eau(:) = fmagic*alb_eau(:) + pmagic
240       alb_eau=MIN(MAX(alb_eau,0.0),1.0)
241!
242    ELSE
243!
244       CALL alboc(REAL(jour),rlat,alb_eau)
245!--ad-hoc correction for model radiative balance tuning
246!--now outside alboc routine
247       alb_eau(:) = fmagic*alb_eau(:) + pmagic
248       alb_eau=MIN(MAX(alb_eau(i),0.04),0.60)
249!
250    ENDIF
251!
252    DO i =1, knon
253      DO  k=1,nsw
254       alb_dir_new(i,k) = alb_eau(knindex(i))
255      ENDDO
256    ENDDO
257!IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions
258!albedo for diffuse radiation is taken the same as for direct radiation
259     alb_dif_new=alb_dir_new
260!IM 09122015 end
261!
262ELSE IF (iflag_albedo==1) THEN
263!--new parametrization of ocean surface albedo by Sunghye Baek
264!--albedo for direct and diffuse radiation are different
265!
266    CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
267!
268!--ad-hoc correction for model radiative balance tuning
269    alb_dir_new(:,:) = fmagic*alb_dir_new(:,:) + pmagic
270    alb_dir_new=MIN(MAX(alb_dir_new,0.0),1.0)
271    alb_dif_new=MIN(MAX(alb_dif_new,0.0),1.0)
272!
273ENDIF
274!albedo SB <<<
275
276!******************************************************************************
277! Calculate the rugosity
278!******************************************************************************
279IF (iflag_z0_oce==0) THEN
280    DO i = 1, knon
281       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
282       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
283            +  0.11*14e-6 / SQRT(cdragm(i) * tmp)
284       z0m(i) = MAX(1.5e-05,z0m(i))
285    ENDDO   
286    z0h(1:knon)=z0m(1:knon) ! En attendant mieux
287
288ELSE IF (iflag_z0_oce==1) THEN
289    DO i = 1, knon
290       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
291       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
292            + 0.11*14e-6 / SQRT(cdragm(i) * tmp)
293       z0m(i) = MAX(1.5e-05,z0m(i))
294       z0h(i)=0.4*14e-6 / SQRT(cdragm(i) * tmp)
295    ENDDO
296ELSE IF (iflag_z0_oce==-1) THEN
297    DO i = 1, knon
298       z0m(i) = z0min
299       z0h(i) = z0min
300    ENDDO
301ELSE
302       CALL abort_physic(modname,'version non prevue',1)
303ENDIF
304!
305!******************************************************************************
306  END SUBROUTINE surf_ocean
307!******************************************************************************
308!
309END MODULE surf_ocean_mod
Note: See TracBrowser for help on using the repository browser.