source: LMDZ5/trunk/libf/phylmd/phyaqua_mod.F90 @ 2225

Last change on this file since 2225 was 2209, checked in by Ehouarn Millour, 10 years ago

Update of the slab ocean by Francis Codron. There are now 3 possibilities for the "version_ocean" slab type:
sicOBS = prescribed ice fraction. Water temperature nearby is set to -1.8°C and cannot become lower.
sicNO = ignore sea ice. One can prescribe a fraction, but the nearby ocean evolves freely, depending on surface fluxes: temperature can go below freezing point or above...
sicINT = interactive sea ice.
EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 29.1 KB
Line 
1MODULE phyaqua_mod
2  ! Routines complementaires pour la physique planetaire.
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE iniaqua(nlon, latfi, lonfi, iflag_phys)
8
9    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10    ! Creation d'un etat initial et de conditions aux limites
11    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
12    ! du modele LMDZ dans sa version terrestre.
13    ! iflag_phys est un parametre qui controle
14    ! iflag_phys = N
15    ! de 100 a 199 : aqua planetes avec SST forcees
16    ! N-100 determine le type de SSTs
17    ! de 200 a 299 : terra planetes avec Ts calcule
18
19    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21    USE comgeomphy, ONLY: rlatd, rlond
22    USE dimphy, ONLY: klon
23    USE surface_data, ONLY: type_ocean, ok_veget
24    USE pbl_surface_mod, ONLY: pbl_surface_init
25    USE fonte_neige_mod, ONLY: fonte_neige_init
26    USE phys_state_var_mod
27    USE control_mod, ONLY: dayref, nday, iphysiq
28    USE indice_sol_mod
29
30    USE ioipsl
31    IMPLICIT NONE
32
33    include "dimensions.h"
34    ! #include "dimphy.h"
35    ! #include "YOMCST.h"
36    include "comconst.h"
37    include "clesphys.h"
38    include "dimsoil.h"
39    include "temps.h"
40
41    INTEGER, INTENT (IN) :: nlon, iflag_phys
42    ! IM ajout latfi, lonfi
43    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
44
45    INTEGER type_profil, type_aqua
46
47    ! Ajouts initialisation des surfaces
48    REAL :: run_off_lic_0(nlon)
49    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
50    REAL :: frugs(nlon, nbsrf)
51    REAL :: agesno(nlon, nbsrf)
52    REAL :: tsoil(nlon, nsoilmx, nbsrf)
53    REAL :: tslab(nlon), seaice(nlon)
54    REAL evap(nlon, nbsrf), fder(nlon)
55
56
57
58    ! Arguments :
59    ! -----------
60
61    ! integer radpas
62    INTEGER it, unit, i, k, itap
63
64    REAL airefi, zcufi, zcvfi
65
66    REAL rugos, albedo
67    REAL tsurf
68    REAL time, timestep, day, day0
69    REAL qsol_f, qsol(nlon)
70    REAL rugsrel(nlon)
71    ! real zmea(nlon),zstd(nlon),zsig(nlon)
72    ! real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
73    ! real rlon(nlon),rlat(nlon)
74    LOGICAL alb_ocean
75    ! integer demih_pas
76
77    CHARACTER *80 ans, file_forctl, file_fordat, file_start
78    CHARACTER *100 file, var
79    CHARACTER *2 cnbl
80
81    REAL phy_nat(nlon, 360)
82    REAL phy_alb(nlon, 360)
83    REAL phy_sst(nlon, 360)
84    REAL phy_bil(nlon, 360)
85    REAL phy_rug(nlon, 360)
86    REAL phy_ice(nlon, 360)
87    REAL phy_fter(nlon, 360)
88    REAL phy_foce(nlon, 360)
89    REAL phy_fsic(nlon, 360)
90    REAL phy_flic(nlon, 360)
91
92    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
93
94    ! intermediate variables to use getin (need to be "save" to be shared by
95    ! all threads)
96    INTEGER, SAVE :: nbapp_rad_omp
97    REAL, SAVE :: co2_ppm_omp, solaire_omp
98    LOGICAL, SAVE :: alb_ocean_omp
99    REAL, SAVE :: rugos_omp
100    ! -------------------------------------------------------------------------
101    ! declaration pour l'appel a phyredem
102    ! -------------------------------------------------------------------------
103
104    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
105    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
106    ! real pbl_tke(nlon,llm,nbsrf)
107    ! real rain_fall(nlon),snow_fall(nlon)
108    ! real solsw(nlon), sollw(nlon),radsol(nlon)
109    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
110    ! real ratqs(nlon,llm)
111    ! real clwcon(nlon,llm)
112
113    INTEGER longcles
114    PARAMETER (longcles=20)
115    REAL clesphy0(longcles)
116
117
118    ! -----------------------------------------------------------------------
119    ! dynamial tendencies :
120    ! ---------------------
121
122    INTEGER l, ierr, aslun
123
124    REAL longitude, latitude
125    REAL paire
126
127    DATA latitude, longitude/48., 0./
128
129    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130    ! INITIALISATIONS
131    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
132
133    ! -----------------------------------------------------------------------
134    ! Initialisations  des constantes
135    ! -------------------------------
136
137
138    type_aqua = iflag_phys/100
139    type_profil = iflag_phys - type_aqua*100
140    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
141
142    IF (klon/=nlon) THEN
143      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
144      STOP 'probleme de dimensions dans iniaqua'
145    END IF
146    CALL phys_state_var_init(read_climoz)
147
148
149    read_climoz = 0
150    day0 = 217.
151    day = day0
152    it = 0
153    time = 0.
154
155    ! IM ajout latfi, lonfi
156    rlatd = latfi
157    rlond = lonfi
158    rlat = rlatd*180./pi
159    rlon = rlond*180./pi
160
161    ! -----------------------------------------------------------------------
162    ! initialisations de la physique
163    ! -----------------------------------------------------------------------
164
165    day_ini = dayref
166    day_end = day_ini + nday
167    airefi = 1.
168    zcufi = 1.
169    zcvfi = 1.
170    !$OMP MASTER
171    nbapp_rad_omp = 24
172    CALL getin('nbapp_rad', nbapp_rad_omp)
173    !$OMP END MASTER
174    !$OMP BARRIER
175    nbapp_rad = nbapp_rad_omp
176
177    ! ---------------------------------------------------------------------
178    ! Creation des conditions aux limites:
179    ! ------------------------------------
180    ! Initialisations des constantes
181    ! Ajouter les manquants dans planete.def... (albedo etc)
182    !$OMP MASTER
183    co2_ppm_omp = 348.
184    CALL getin('co2_ppm', co2_ppm_omp)
185    solaire_omp = 1365.
186    CALL getin('solaire', solaire_omp)
187    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
188    ! type_aqua
189    alb_ocean_omp = .TRUE.
190    CALL getin('alb_ocean', alb_ocean_omp)
191    !$OMP END MASTER
192    !$OMP BARRIER
193    co2_ppm = co2_ppm_omp
194    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
195    solaire = solaire_omp
196    WRITE (*, *) 'iniaqua: solaire=', solaire
197    alb_ocean = alb_ocean_omp
198    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
199
200    radsol = 0.
201    qsol_f = 10.
202
203    ! Conditions aux limites:
204    ! -----------------------
205
206    qsol(:) = qsol_f
207    rugsrel = 0.0 ! (rugsrel = rugoro)
208    rugoro = 0.0
209    u_ancien = 0.0
210    v_ancien = 0.0
211    agesno = 50.0
212    ! Relief plat
213    zmea = 0.
214    zstd = 0.
215    zsig = 0.
216    zgam = 0.
217    zthe = 0.
218    zpic = 0.
219    zval = 0.
220
221    ! Une seule surface
222    pctsrf = 0.
223    IF (type_aqua==1) THEN
224      rugos = 1.E-4
225      albedo = 0.19
226      pctsrf(:, is_oce) = 1.
227    ELSE IF (type_aqua==2) THEN
228      rugos = 0.03
229      albedo = 0.1
230      pctsrf(:, is_ter) = 1.
231    END IF
232
233    !$OMP MASTER
234    rugos_omp = rugos
235    CALL getin('rugos', rugos_omp)
236    !$OMP END MASTER
237    !$OMP BARRIER
238    rugos = rugos_omp
239    WRITE (*, *) 'iniaqua: rugos=', rugos
240    zmasq(:) = pctsrf(:, is_ter)
241
242    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
243    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
244
245    ! Si alb_ocean on calcule un albedo oceanique moyen
246    ! if (alb_ocean) then
247    ! Voir pourquoi on avait ca.
248    ! CALL ini_alb_oce(phy_alb)
249    ! else
250    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
251    ! endif !alb_ocean
252
253    DO i = 1, 360
254      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
255      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
256
257      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
258      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
259      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
260      phy_ice(:, i) = 0.0 ! fraction de glace (?)
261      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
262      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
263      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
264      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
265    END DO
266    ! IM calcul profil sst
267    CALL profil_sst(nlon, rlatd, type_profil, phy_sst)
268
269    CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
270      phy_fter, phy_foce, phy_flic, phy_fsic)
271
272
273    ! ---------------------------------------------------------------------
274    ! Ecriture de l'etat initial:
275    ! ---------------------------
276
277
278    ! Ecriture etat initial physique
279
280    timestep = dtvr*float(iphysiq)
281    radpas = nint(daysec/timestep/float(nbapp_rad))
282
283    DO i = 1, longcles
284      clesphy0(i) = 0.
285    END DO
286    clesphy0(1) = float(iflag_con)
287    clesphy0(2) = float(nbapp_rad)
288    ! IF( cycle_diurne  ) clesphy0(3) =  1.
289    clesphy0(3) = 1. ! cycle_diurne
290    clesphy0(4) = 1. ! soil_model
291    clesphy0(5) = 1. ! new_oliq
292    clesphy0(6) = 0. ! ok_orodr
293    clesphy0(7) = 0. ! ok_orolf
294    clesphy0(8) = 0. ! ok_limitvrai
295
296
297    ! =======================================================================
298    ! Profils initiaux
299    ! =======================================================================
300
301    ! On initialise les temperatures de surfaces comme les sst
302    DO i = 1, nlon
303      ftsol(i, :) = phy_sst(i, 1)
304      tsoil(i, :, :) = phy_sst(i, 1)
305      tslab(i) = phy_sst(i, 1)
306    END DO
307
308    falbe(:, :) = albedo
309    falblw(:, :) = albedo
310    rain_fall(:) = 0.
311    snow_fall(:) = 0.
312    solsw(:) = 0.
313    sollw(:) = 0.
314    radsol(:) = 0.
315
316    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
317    ! intialisation bidon mais pas grave
318    t_ancien(:, :) = 0.
319    q_ancien(:, :) = 0.
320    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
321    rnebcon = 0.
322    ratqs = 0.
323    clwcon = 0.
324    pbl_tke = 1.E-8
325
326    ! variables supplementaires pour appel a plb_surface_init
327    fder(:) = 0.
328    seaice(:) = 0.
329    run_off_lic_0 = 0.
330    evap = 0.
331
332
333    ! Initialisations necessaires avant phyredem
334    type_ocean = 'force'
335    CALL fonte_neige_init(run_off_lic_0)
336    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
337    snsrf(:, :) = 0. ! couverture de neige des sous surface
338    frugs(:, :) = rugos ! couverture de neige des sous surface
339
340
341    CALL pbl_surface_init(qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, &
342      tsoil)
343
344    PRINT *, 'iniaqua: before phyredem'
345
346    falb1 = albedo
347    falb2 = albedo
348    zmax0 = 0.
349    f0 = 0.
350    sig1 = 0.
351    w01 = 0.
352    wake_deltat = 0.
353    wake_deltaq = 0.
354    wake_s = 0.
355    wake_cstar = 0.
356    wake_pe = 0.
357    wake_fip = 0.
358    fm_therm = 0.
359    entr_therm = 0.
360    detr_therm = 0.
361
362
363    CALL phyredem('startphy.nc')
364
365    PRINT *, 'iniaqua: after phyredem'
366    CALL phys_state_var_end
367
368    RETURN
369  END SUBROUTINE iniaqua
370
371
372  ! ====================================================================
373  ! ====================================================================
374  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
375    USE dimphy
376    IMPLICIT NONE
377    ! ====================================================================
378    ! =============================================================
379    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
380    ! Auteur : A. Campoy et F. Hourdin
381    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
382    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
383    ! connaissant la declinaison, la latitude et la longitude.
384
385    ! Dans cette version particuliere, on calcule le rayonnement
386    ! moyen sur l'année à chaque latitude.
387    ! angle zenithal calculé pour obtenir un
388    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
389    ! en moyenne annuelle.
390    ! Spécifique de la terre. Utilisé pour les aqua planetes.
391
392    ! Rque   : Different de la routine angle en ce sens que zenang
393    ! fournit des moyennes de pmu0 et non des valeurs
394    ! instantanees, du coup frac prend toutes les valeurs
395    ! entre 0 et 1.
396    ! Date   : premiere version le 13 decembre 1994
397    ! revu pour  GCM  le 30 septembre 1996
398    ! ===============================================================
399    ! longi----INPUT : la longitude vraie de la terre dans son plan
400    ! solaire a partir de l'equinoxe de printemps (degre)
401    ! gmtime---INPUT : temps universel en fraction de jour
402    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
403    ! lat------INPUT : latitude en degres
404    ! long-----INPUT : longitude en degres
405    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
406    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
407    ! ================================================================
408    include "YOMCST.h"
409    ! ================================================================
410    LOGICAL cycle_diurne
411    REAL gmtime
412    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
413    ! ================================================================
414    INTEGER i
415    REAL gmtime1, gmtime2
416    REAL pi_local
417
418
419    REAL rmu0m(klon), rmu0a(klon)
420
421
422    pi_local = 4.0*atan(1.0)
423
424    ! ================================================================
425    ! Calcul de l'angle zenithal moyen sur la journee
426    ! ================================================================
427
428    DO i = 1, klon
429      fract(i) = 1.
430      ! Calcule du flux moyen
431      IF (abs(rlat(i))<=28.75) THEN
432        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
433      ELSE IF (abs(rlat(i))<=43.75) THEN
434        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
435      ELSE IF (abs(rlat(i))<=71.25) THEN
436        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
437      ELSE
438        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
439      END IF
440    END DO
441
442    ! ================================================================
443    ! Avec ou sans cycle diurne
444    ! ================================================================
445
446    IF (cycle_diurne) THEN
447
448      ! On redecompose flux  au sommet suivant un cycle diurne idealise
449      ! identique a toutes les latitudes.
450
451      DO i = 1, klon
452        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
453        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
454          rmu0a(i)/sqrt(2.)
455      END DO
456
457      DO i = 1, klon
458        IF (rmu0(i)<=0.) THEN
459          rmu0(i) = 0.
460          fract(i) = 0.
461        ELSE
462          fract(i) = 1.
463        END IF
464      END DO
465
466      ! Affichage de l'angel zenitale
467      ! print*,'************************************'
468      ! print*,'************************************'
469      ! print*,'************************************'
470      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
471      ! print*,'rmu0m=',rmu0m(i)
472      ! print*,'rmu0a=',rmu0a(i)
473      ! print*,'rmu0=',rmu0(i)
474
475    ELSE
476
477      DO i = 1, klon
478        fract(i) = 0.5
479        rmu0(i) = rmu0m(i)*2.
480      END DO
481
482    END IF
483
484    RETURN
485  END SUBROUTINE zenang_an
486
487  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488
489  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
490      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
491
492    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_omp_root
493    USE mod_grid_phy_lmdz, ONLY: klon_glo
494    USE mod_phys_lmdz_transfert_para, ONLY: gather
495    IMPLICIT NONE
496    ! #include "dimensions.h"
497    ! #include "dimphy.h"
498    include "netcdf.inc"
499
500    INTEGER, INTENT (IN) :: klon
501    REAL, INTENT (IN) :: phy_nat(klon, 360)
502    REAL, INTENT (IN) :: phy_alb(klon, 360)
503    REAL, INTENT (IN) :: phy_sst(klon, 360)
504    REAL, INTENT (IN) :: phy_bil(klon, 360)
505    REAL, INTENT (IN) :: phy_rug(klon, 360)
506    REAL, INTENT (IN) :: phy_ice(klon, 360)
507    REAL, INTENT (IN) :: phy_fter(klon, 360)
508    REAL, INTENT (IN) :: phy_foce(klon, 360)
509    REAL, INTENT (IN) :: phy_flic(klon, 360)
510    REAL, INTENT (IN) :: phy_fsic(klon, 360)
511
512    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
513      ! on the whole physics grid
514    INTEGER :: k
515    INTEGER ierr
516    INTEGER dimfirst(3)
517    INTEGER dimlast(3)
518
519    INTEGER nid, ndim, ntim
520    INTEGER dims(2), debut(2), epais(2)
521    INTEGER id_tim
522    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
523    INTEGER id_fter, id_foce, id_fsic, id_flic
524
525    IF (is_mpi_root .AND. is_omp_root) THEN
526
527      PRINT *, 'writelim: Ecriture du fichier limit'
528
529      ierr = nf_create('limit.nc', nf_clobber, nid)
530
531      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
532        'Fichier conditions aux limites')
533      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
534      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
535      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
536
537      dims(1) = ndim
538      dims(2) = ntim
539
540#ifdef NC_DOUBLE
541      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
542#else
543      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
544#endif
545      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
546
547#ifdef NC_DOUBLE
548      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
549#else
550      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
551#endif
552      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
553        'Nature du sol (0,1,2,3)')
554
555#ifdef NC_DOUBLE
556      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
557#else
558      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
559#endif
560      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
561        'Temperature superficielle de la mer')
562
563#ifdef NC_DOUBLE
564      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
565#else
566      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
567#endif
568      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
569        'Reference flux de chaleur au sol')
570
571#ifdef NC_DOUBLE
572      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
573#else
574      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
575#endif
576      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
577
578#ifdef NC_DOUBLE
579      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
580#else
581      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
582#endif
583      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
584
585#ifdef NC_DOUBLE
586      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
587#else
588      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
589#endif
590      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
591#ifdef NC_DOUBLE
592      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
593#else
594      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
595#endif
596      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
597#ifdef NC_DOUBLE
598      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
599#else
600      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
601#endif
602      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
603#ifdef NC_DOUBLE
604      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
605#else
606      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
607#endif
608      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
609
610      ierr = nf_enddef(nid)
611      IF (ierr/=nf_noerr) THEN
612        WRITE (*, *) 'writelim error: failed to end define mode'
613        WRITE (*, *) nf_strerror(ierr)
614      END IF
615
616
617      ! write the 'times'
618      DO k = 1, 360
619#ifdef NC_DOUBLE
620        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
621#else
622        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
623#endif
624        IF (ierr/=nf_noerr) THEN
625          WRITE (*, *) 'writelim error with temps(k),k=', k
626          WRITE (*, *) nf_strerror(ierr)
627        END IF
628      END DO
629
630    END IF ! of if (is_mpi_root.and.is_omp_root)
631
632    ! write the fields, after having collected them on master
633
634    CALL gather(phy_nat, phy_glo)
635    IF (is_mpi_root .AND. is_omp_root) THEN
636#ifdef NC_DOUBLE
637      ierr = nf_put_var_double(nid, id_nat, phy_glo)
638#else
639      ierr = nf_put_var_real(nid, id_nat, phy_glo)
640#endif
641      IF (ierr/=nf_noerr) THEN
642        WRITE (*, *) 'writelim error with phy_nat'
643        WRITE (*, *) nf_strerror(ierr)
644      END IF
645    END IF
646
647    CALL gather(phy_sst, phy_glo)
648    IF (is_mpi_root .AND. is_omp_root) THEN
649#ifdef NC_DOUBLE
650      ierr = nf_put_var_double(nid, id_sst, phy_glo)
651#else
652      ierr = nf_put_var_real(nid, id_sst, phy_glo)
653#endif
654      IF (ierr/=nf_noerr) THEN
655        WRITE (*, *) 'writelim error with phy_sst'
656        WRITE (*, *) nf_strerror(ierr)
657      END IF
658    END IF
659
660    CALL gather(phy_bil, phy_glo)
661    IF (is_mpi_root .AND. is_omp_root) THEN
662#ifdef NC_DOUBLE
663      ierr = nf_put_var_double(nid, id_bils, phy_glo)
664#else
665      ierr = nf_put_var_real(nid, id_bils, phy_glo)
666#endif
667      IF (ierr/=nf_noerr) THEN
668        WRITE (*, *) 'writelim error with phy_bil'
669        WRITE (*, *) nf_strerror(ierr)
670      END IF
671    END IF
672
673    CALL gather(phy_alb, phy_glo)
674    IF (is_mpi_root .AND. is_omp_root) THEN
675#ifdef NC_DOUBLE
676      ierr = nf_put_var_double(nid, id_alb, phy_glo)
677#else
678      ierr = nf_put_var_real(nid, id_alb, phy_glo)
679#endif
680      IF (ierr/=nf_noerr) THEN
681        WRITE (*, *) 'writelim error with phy_alb'
682        WRITE (*, *) nf_strerror(ierr)
683      END IF
684    END IF
685
686    CALL gather(phy_rug, phy_glo)
687    IF (is_mpi_root .AND. is_omp_root) THEN
688#ifdef NC_DOUBLE
689      ierr = nf_put_var_double(nid, id_rug, phy_glo)
690#else
691      ierr = nf_put_var_real(nid, id_rug, phy_glo)
692#endif
693      IF (ierr/=nf_noerr) THEN
694        WRITE (*, *) 'writelim error with phy_rug'
695        WRITE (*, *) nf_strerror(ierr)
696      END IF
697    END IF
698
699    CALL gather(phy_fter, phy_glo)
700    IF (is_mpi_root .AND. is_omp_root) THEN
701#ifdef NC_DOUBLE
702      ierr = nf_put_var_double(nid, id_fter, phy_glo)
703#else
704      ierr = nf_put_var_real(nid, id_fter, phy_glo)
705#endif
706      IF (ierr/=nf_noerr) THEN
707        WRITE (*, *) 'writelim error with phy_fter'
708        WRITE (*, *) nf_strerror(ierr)
709      END IF
710    END IF
711
712    CALL gather(phy_foce, phy_glo)
713    IF (is_mpi_root .AND. is_omp_root) THEN
714#ifdef NC_DOUBLE
715      ierr = nf_put_var_double(nid, id_foce, phy_glo)
716#else
717      ierr = nf_put_var_real(nid, id_foce, phy_glo)
718#endif
719      IF (ierr/=nf_noerr) THEN
720        WRITE (*, *) 'writelim error with phy_foce'
721        WRITE (*, *) nf_strerror(ierr)
722      END IF
723    END IF
724
725    CALL gather(phy_fsic, phy_glo)
726    IF (is_mpi_root .AND. is_omp_root) THEN
727#ifdef NC_DOUBLE
728      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
729#else
730      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
731#endif
732      IF (ierr/=nf_noerr) THEN
733        WRITE (*, *) 'writelim error with phy_fsic'
734        WRITE (*, *) nf_strerror(ierr)
735      END IF
736    END IF
737
738    CALL gather(phy_flic, phy_glo)
739    IF (is_mpi_root .AND. is_omp_root) THEN
740#ifdef NC_DOUBLE
741      ierr = nf_put_var_double(nid, id_flic, phy_glo)
742#else
743      ierr = nf_put_var_real(nid, id_flic, phy_glo)
744#endif
745      IF (ierr/=nf_noerr) THEN
746        WRITE (*, *) 'writelim error with phy_flic'
747        WRITE (*, *) nf_strerror(ierr)
748      END IF
749    END IF
750
751    ! close file:
752    IF (is_mpi_root .AND. is_omp_root) THEN
753      ierr = nf_close(nid)
754    END IF
755
756  END SUBROUTINE writelim
757
758  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
759
760  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
761    USE dimphy
762    IMPLICIT NONE
763
764    INTEGER nlon, type_profil, i, k, j
765    REAL :: rlatd(nlon), phy_sst(nlon, 360)
766    INTEGER imn, imx, amn, amx, kmn, kmx
767    INTEGER p, pplus, nlat_max
768    PARAMETER (nlat_max=72)
769    REAL x_anom_sst(nlat_max)
770
771    IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua'
772    WRITE (*, *) ' profil_sst: type_profil=', type_profil
773    DO i = 1, 360
774      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
775
776      ! Rajout fbrlmd
777
778      IF (type_profil==1) THEN
779        ! Méthode 1 "Control" faible plateau à l'Equateur
780        DO j = 1, klon
781          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
782          ! PI/3=1.047197551
783          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
784            phy_sst(j, i) = 273.
785          END IF
786        END DO
787      END IF
788      IF (type_profil==2) THEN
789        ! Méthode 2 "Flat" fort plateau à l'Equateur
790        DO j = 1, klon
791          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
792          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
793            phy_sst(j, i) = 273.
794          END IF
795        END DO
796      END IF
797
798
799      IF (type_profil==3) THEN
800        ! Méthode 3 "Qobs" plateau réel à l'Equateur
801        DO j = 1, klon
802          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
803            rlatd(j))**4)
804          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
805            phy_sst(j, i) = 273.
806          END IF
807        END DO
808      END IF
809
810      IF (type_profil==4) THEN
811        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
812        DO j = 1, klon
813          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
814            rlatd(j))**4)
815          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
816            phy_sst(j, i) = 273.
817          END IF
818        END DO
819      END IF
820
821      IF (type_profil==5) THEN
822        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
823        DO j = 1, klon
824          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
825            *rlatd(j))**4)
826          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
827            phy_sst(j, i) = 273. + 2.
828          END IF
829
830        END DO
831      END IF
832
833      IF (type_profil==6) THEN
834        ! Méthode 6 "cst" valeur constante de SST
835        DO j = 1, klon
836          phy_sst(j, i) = 288.
837        END DO
838      END IF
839
840
841      IF (type_profil==7) THEN
842        ! Méthode 7 "cst" valeur constante de SST +2
843        DO j = 1, klon
844          phy_sst(j, i) = 288. + 2.
845        END DO
846      END IF
847
848      p = 0
849      IF (type_profil==8) THEN
850        ! Méthode 8 profil anomalies SST du modèle couplé AR4
851        DO j = 1, klon
852          IF (rlatd(j)==rlatd(j-1)) THEN
853            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
854              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
855            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
856              phy_sst(j, i) = 273. + x_anom_sst(pplus)
857            END IF
858          ELSE
859            p = p + 1
860            pplus = 73 - p
861            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
862              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
863            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
864              phy_sst(j, i) = 273. + x_anom_sst(pplus)
865            END IF
866            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
867          END IF
868        END DO
869      END IF
870
871      IF (type_profil==9) THEN
872        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
873        DO j = 1, klon
874          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
875            *rlatd(j))**4)
876          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
877            phy_sst(j, i) = 273. - 2.
878          END IF
879        END DO
880      END IF
881
882
883      IF (type_profil==10) THEN
884        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
885        DO j = 1, klon
886          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
887            *rlatd(j))**4)
888          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
889            phy_sst(j, i) = 273. + 4.
890          END IF
891        END DO
892      END IF
893
894      IF (type_profil==11) THEN
895        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
896        DO j = 1, klon
897          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
898            rlatd(j))**4)
899          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
900            phy_sst(j, i) = 273.
901          END IF
902        END DO
903      END IF
904
905      IF (type_profil==12) THEN
906        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
907        DO j = 1, klon
908          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
909            *rlatd(j))**4)
910          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
911            phy_sst(j, i) = 273. + 4.
912          END IF
913        END DO
914      END IF
915
916      IF (type_profil==13) THEN
917        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
918        DO j = 1, klon
919          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
920            rlatd(j))**4)
921          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
922            phy_sst(j, i) = 273.
923          END IF
924        END DO
925      END IF
926
927      IF (type_profil==14) THEN
928        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
929        DO j = 1, klon
930          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
931            *rlatd(j))**4)
932          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
933            phy_sst(j, i) = 273.
934          END IF
935        END DO
936      END IF
937
938      if (type_profil.EQ.20) then
939      print*,'Profile SST 20'
940!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
941
942      do j=1,klon
943        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
944      enddo
945      endif
946
947      if (type_profil.EQ.21) then
948      print*,'Profile SST 21'
949!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
950      do j=1,klon
951        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
952      enddo
953      endif
954
955
956
957    END DO
958
959    ! IM beg : verif profil SST: phy_sst
960    amn = min(phy_sst(1,1), 1000.)
961    amx = max(phy_sst(1,1), -1000.)
962    imn = 1
963    kmn = 1
964    imx = 1
965    kmx = 1
966    DO k = 1, 360
967      DO i = 2, nlon
968        IF (phy_sst(i,k)<amn) THEN
969          amn = phy_sst(i, k)
970          imn = i
971          kmn = k
972        END IF
973        IF (phy_sst(i,k)>amx) THEN
974          amx = phy_sst(i, k)
975          imx = i
976          kmx = k
977        END IF
978      END DO
979    END DO
980
981    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
982    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
983    ! IM end : verif profil SST: phy_sst
984
985    RETURN
986  END SUBROUTINE profil_sst
987
988END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.