source: trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 @ 3322

Last change on this file since 3322 was 3100, checked in by bhatnags, 20 months ago

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM
File size: 13.8 KB
Line 
1module phyetat0_mod
2
3implicit none
4
5contains
6
7subroutine phyetat0 (startphy_file, &
8                     ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
9                     day_ini,time,tsurf,tsoil, &
10                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
11                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
12
13! to use  'getin_p'
14      use ioipsl_getin_p_mod, only: getin_p
15!!
16  use write_field_phy, only: Writefield_phy
17!!
18  use tabfi_mod, only: tabfi
19  USE tracer_h, ONLY: noms, igcm_h2o_vap
20  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
21  use iostart, only: nid_start, open_startphy, close_startphy, &
22                     get_field, get_var, inquire_field, &
23                     inquire_dimension, inquire_dimension_length
24!  use slab_ice_h, only: noceanmx
25  USE ocean_slab_mod, ONLY: nslay
26  use callkeys_mod, only: CLFvarying,surfalbedo,surfemis, callsoil
27  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
28                                east_gwstress, west_gwstress
29  implicit none
30
31!======================================================================
32!  Arguments:
33!  ---------
34!  inputs:
35  logical,intent(in) :: startphy_file ! .true. if reading start file
36  integer,intent(in) :: ngrid
37  integer,intent(in) :: nlayer
38  character*(*),intent(in) :: fichnom ! "startfi.nc" file
39  integer,intent(in) :: tab0
40  integer,intent(in) :: Lmodif
41  integer,intent(in) :: nsoil ! # of soil layers
42  integer,intent(in) :: nq
43  integer,intent(out) :: day_ini
44  real,intent(out) :: time
45
46!  outputs:
47  real,intent(out) :: tsurf(ngrid) ! surface temperature
48  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
49  real,intent(out) :: emis(ngrid) ! surface emissivity
50  real,intent(out) :: q2(ngrid,nlayer+1) !
51  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
52! real co2ice(ngrid) ! co2 ice cover
53  real,intent(out) :: cloudfrac(ngrid,nlayer)
54  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
55  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,nslay) 
56  real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
57  real,intent(out) :: rnat(ngrid)
58!  real,intent(out) :: tice(ngrid)
59
60!======================================================================
61!  Local variables:
62
63!      INTEGER radpas
64!      REAL co2_ppm
65!      REAL solaire
66
67      real xmin,xmax ! to display min and max of a field
68!
69      INTEGER ig,iq,lmax
70      INTEGER nid, nvarid
71      INTEGER ierr, i, nsrf
72!      integer isoil
73!      INTEGER length
74!      PARAMETER (length=100)
75      CHARACTER*7 str7
76      CHARACTER*2 str2
77      CHARACTER*1 yes
78!
79      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
80      INTEGER nqold
81
82! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
83!      logical :: oldtracernames=.false.
84      integer :: count
85      character(len=30) :: txt ! to store some text
86     
87      INTEGER :: indextime=1 ! index of selected time, default value=1
88      logical :: found
89     
90      character(len=8) :: modname="phyetat0"
91
92!
93! ALLOCATE ARRAYS IN surfdat_h
94!
95IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
96IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
97IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
98IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
99IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
100IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
101IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
102
103if (startphy_file) then
104  ! open physics initial state file:
105  call open_startphy(fichnom)
106
107  ! possibility to modify tab_cntrl in tabfi
108  write(*,*)
109  write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
110  call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
111                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
112
113else ! "academic" initialization of planetary parameters via tabfi
114  call tabfi (ngrid,0,0,0,day_ini,lmax,p_rad, &
115                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
116endif ! of if (startphy_file)
117
118if (startphy_file) then
119  ! Load surface geopotential:
120  call get_field("phisfi",phisfi,found)
121  if (.not.found) then
122    call abort_physic(modname,"Failed loading <phisfi>",1)
123  endif
124else
125  phisfi(:)=0.
126endif ! of if (startphy_file)
127write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
128               minval(phisfi), maxval(phisfi)
129
130if (startphy_file) then
131  ! Load bare ground albedo:
132  call get_field("albedodat",albedodat,found)
133  if (.not.found) then
134    call abort_physic(modname,"Failed loading <albedodat>",1)
135  endif
136else
137  ! If no startfi file, use parameter surfalbedo in def file
138  surfalbedo=0.5
139  call getin_p("surfalbedo",surfalbedo)
140  print*,"surfalbedo",surfalbedo
141  albedodat(:)=surfalbedo
142endif ! of if (startphy_file)
143write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
144             minval(albedodat), maxval(albedodat)
145
146! ZMEA
147if (startphy_file) then
148  call get_field("ZMEA",zmea,found)
149  if (.not.found) then
150    call abort_physic(modname,"Failed loading <ZMEA>",1)
151  endif
152else
153  zmea(:)=0.
154endif ! of if (startphy_file)
155write(*,*) "phyetat0: <ZMEA> range:", &
156             minval(zmea), maxval(zmea)
157
158! ZSTD
159if (startphy_file) then
160  call get_field("ZSTD",zstd,found)
161  if (.not.found) then
162    call abort_physic(modname,"Failed loading <ZSTD>",1)
163  endif
164else
165  zstd(:)=0.
166endif ! of if (startphy_file)
167write(*,*) "phyetat0: <ZSTD> range:", &
168             minval(zstd), maxval(zstd)
169
170! ZSIG
171if (startphy_file) then
172  call get_field("ZSIG",zsig,found)
173  if (.not.found) then
174    call abort_physic(modname,"Failed loading <ZSIG>",1)
175  endif
176else
177  zsig(:)=0.
178endif ! of if (startphy_file)
179write(*,*) "phyetat0: <ZSIG> range:", &
180             minval(zsig), maxval(zsig)
181
182! ZGAM
183if (startphy_file) then
184  call get_field("ZGAM",zgam,found)
185  if (.not.found) then
186    call abort_physic(modname,"Failed loading <ZGAM>",1)
187  endif
188else
189  zgam(:)=0.
190endif ! of if (startphy_file)
191write(*,*) "phyetat0: <ZGAM> range:", &
192             minval(zgam), maxval(zgam)
193
194! ZTHE
195if (startphy_file) then
196  call get_field("ZTHE",zthe,found)
197  if (.not.found) then
198    call abort_physic(modname,"Failed loading <ZTHE>",1)
199  endif
200else
201  zthe(:)=0.
202endif ! of if (startphy_file)
203write(*,*) "phyetat0: <ZTHE> range:", &
204             minval(zthe), maxval(zthe)
205
206! Surface temperature :
207if (startphy_file) then
208  call get_field("tsurf",tsurf,found,indextime)
209  if (.not.found) then
210    call abort_physic(modname,"Failed loading <tsurf>",1)
211  endif
212else
213  tsurf(:)=0. ! will be updated afterwards in physiq !
214endif ! of if (startphy_file)
215write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
216             minval(tsurf), maxval(tsurf)
217
218! Surface emissivity
219if (startphy_file) then
220  call get_field("emis",emis,found,indextime)
221  if (.not.found) then
222    call abort_physic(modname,"Failed loading <emis>",1)
223  endif
224else
225  ! If no startfi file, use parameter surfemis in def file
226  surfemis=1.0
227  call getin_p("surfemis",surfemis)
228  print*,"surfemis",surfemis
229  emis(:)=surfemis
230endif ! of if (startphy_file)
231write(*,*) "phyetat0: Surface emissivity <emis> range:", &
232             minval(emis), maxval(emis)
233
234! Cloud fraction (added by BC 2010)
235if (CLFvarying) then
236  if (startphy_file) then
237    call get_field("cloudfrac",cloudfrac,found,indextime)
238    if (.not.found) then
239      call abort_physic(modname,"Failed loading <cloudfrac>",1)
240    endif
241  else
242    cloudfrac(:,:)=0.0
243  endif ! of if (startphy_file)
244  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
245               minval(cloudfrac), maxval(cloudfrac)
246else
247  cloudfrac(:,:)=0.0
248endif ! of if (CLFvarying)
249
250! Total cloud fraction (added by BC 2010)
251if (CLFvarying) then
252  if (startphy_file) then
253    call get_field("totcloudfrac",totcloudfrac,found,indextime)
254    if (.not.found) then
255      call abort_physic(modname,"Failed loading <totcloudfrac>",1)
256    endif
257  else
258    totcloudfrac(:)=0.0
259  endif ! of if (startphy_file)
260  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
261             minval(totcloudfrac), maxval(totcloudfrac)
262else
263  totcloudfrac(:)=0.0
264endif ! of if (CLFvarying)
265
266! Height of oceanic ice (added by BC 2010)
267if (startphy_file) then
268  call get_field("hice",hice,found,indextime)
269  if (.not.found) then
270    write(*,*) "phyetat0: Failed loading <hice>"
271    !  call abort
272    hice(:)=0.
273  endif
274else
275  hice(:)=0.
276endif ! of if (startphy_file)
277write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
278             minval(hice), maxval(hice)
279
280! SLAB OCEAN (added by BC 2014)
281if (startphy_file) then
282  ! nature of the surface
283  call get_field("rnat",rnat,found,indextime)
284  if (.not.found) then
285    write(*,*) "phyetat0: Failed loading <rnat>"
286        rnat(1:ngrid)=1.
287  else
288      do ig=1,ngrid
289        if((nint(rnat(ig)).eq.2).or.(nint(rnat(ig)).eq.0))then
290          rnat(ig)=0.
291        else
292          rnat(ig)=1.
293        endif     
294      enddo
295  endif ! of if (.not.found)
296else
297  rnat(:)=1.
298endif ! of if (startphy_file)
299write(*,*) "phyetat0: Nature of surface <rnat> range:", &
300             minval(rnat), maxval(rnat)
301
302if (startphy_file) then
303  ! Pourcentage of sea ice cover
304  call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
305  if (.not.found) then
306    write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
307    pctsrf_sic(1:ngrid)=0.
308  endif
309else
310  pctsrf_sic(:)=0.
311endif ! of if (startphy_file)
312write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
313             minval(pctsrf_sic), maxval(pctsrf_sic)
314
315if (startphy_file) then
316  ! Slab ocean temperature (2 layers)
317  call get_field("tslab",tslab,found,indextime)
318  if (.not.found) then
319    write(*,*) "phyetat0: Failed loading <tslab>"
320    do iq=1,nslay
321      tslab(1:ngrid,iq)=tsurf(1:ngrid)
322    enddo
323  endif
324else
325  do iq=1,nslay
326    tslab(1:ngrid,iq)=tsurf(1:ngrid)
327  enddo
328endif ! of if (startphy_file)
329write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
330             minval(tslab), maxval(tslab)
331
332if (startphy_file) then
333  ! Oceanic ice temperature
334  call get_field("tsea_ice",tsea_ice,found,indextime)
335  if (.not.found) then
336    write(*,*) "phyetat0: Failed loading <tsea_ice>"
337    tsea_ice(1:ngrid)=273.15-1.8
338  endif
339else
340  tsea_ice(1:ngrid)=273.15-1.8
341endif ! of if (startphy_file)
342write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
343             minval(tsea_ice), maxval(tsea_ice)
344
345!if (startphy_file) then
346  ! Oceanic ice temperature
347!  call get_field("tice",tice,found,indextime)
348!  if (.not.found) then
349!    write(*,*) "phyetat0: Failed loading <tice>"
350!    tice(1:ngrid)=273.15-1.8
351!  endif
352!else
353!  tice(1:ngrid)=273.15-1.8
354!endif ! of if (startphy_file)
355!write(*,*) "phyetat0: Oceanic ice temperature <tice> range:", &
356!             minval(tice), maxval(tice)
357
358if (startphy_file) then
359  !  Oceanic ice quantity (kg/m^2)
360  call get_field("sea_ice",sea_ice,found,indextime)
361  if (.not.found) then
362    write(*,*) "phyetat0: Failed loading <sea_ice>"
363    sea_ice(1:ngrid)=0.
364  endif
365else
366  sea_ice(1:ngrid)=0.
367endif ! of if (startphy_file)
368write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
369             minval(sea_ice), maxval(sea_ice)
370
371
372! pbl wind variance
373if (startphy_file) then
374  call get_field("q2",q2,found,indextime)
375  if (.not.found) then
376    call abort_physic(modname,"Failed loading <q2>",1)
377  endif
378else
379  q2(:,:)=0.
380endif ! of if (startphy_file)
381write(*,*) "phyetat0: PBL wind variance <q2> range:", &
382             minval(q2), maxval(q2)
383
384! tracer on surface
385if (nq.ge.1) then
386  do iq=1,nq
387    txt=noms(iq)
388    if (startphy_file) then
389      call get_field(txt,qsurf(:,iq),found,indextime)
390      if (.not.found) then
391        write(*,*) "phyetat0: Failed loading <",trim(txt),">"
392        write(*,*) "         ",trim(txt)," is set to zero"
393        qsurf(:,iq) = 0.
394      endif
395    else
396      qsurf(:,iq)=0.
397    endif ! of if (startphy_file)
398    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
399                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
400  enddo! of do iq=1,nq
401endif ! of if (nq.ge.1)
402
403!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
404
405if (startphy_file) then
406  ! Call to soil_settings, in order to read soil temperatures,
407  ! as well as thermal inertia and volumetric heat capacity
408  if (callsoil) then
409    call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
410  endif
411endif ! of if (startphy_file)
412
413! Non-orographic gravity waves
414if (startphy_file) then
415   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
416   if (.not.found) then
417      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
418      du_nonoro_gwd(:,:)=0.
419   else
420      write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
421      write(*,*) " <du_nonoro_gwd> range:", &
422             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
423   endif
424endif ! of if (startphy_file)
425if (startphy_file) then
426   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
427   if (.not.found) then
428      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
429      dv_nonoro_gwd(:,:)=0.
430   else
431      write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
432      write(*,*) " <dv_nonoro_gwd> range:", &
433             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
434   endif
435endif ! of if (startphy_file)
436if (startphy_file) then
437   call get_field("east_gwstress",east_gwstress,found,indextime)
438   if (.not.found) then
439      write(*,*) "phyetat0: <east_gwstress> not in file"
440      east_gwstress(:,:)=0.
441   else
442      write(*,*) "phyetat0: Memory of eastward stress profile due to non-orographic GW"
443      write(*,*) " <east_gwstress> range:", &
444             minval(east_gwstress), maxval(east_gwstress)
445   endif
446endif ! of if (startphy_file)
447if (startphy_file) then
448   call get_field("west_gwstress",west_gwstress,found,indextime)
449   if (.not.found) then
450      write(*,*) "phyetat0: <west_gwstress> not in file"
451      west_gwstress(:,:)=0.
452   else
453      write(*,*) "phyetat0: Memory of westward stress profile due to non-orographic GW"
454      write(*,*) " <west_gwstress> range:", &
455             minval(west_gwstress), maxval(west_gwstress)
456   endif
457endif ! of if (startphy_file)
458!
459! close file:
460!
461if (startphy_file) call close_startphy
462
463end subroutine phyetat0
464
465end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.